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PROJECT SUMMARY 


At the start of our work in June 1971, three 
principal tasks were assigned. Slightly restated, 
these tasks were: 

(i) Develop a deterministic algorithm for 
guidance to rendezvous with comets and asteroids 
that can handle expected large target ephemeris 
errors . 

(ii) Define the problem of determination of 
rotational state of a tumbling asteroid or cometary 
nucleus and develop possible schemes for this 
determination . 

(iii) Investigate possible contact rendezvous 
schemes including the "harpoon" technique. 

During the first year, a detailed investigation 
of a rendezvous guidance technique based on 
"encounter theory" was conducted. The definition 
and formulation of the tumbling problem was made 
and several possible algorithms phrased. A first 
investigation of the harpoon problem was conducted 
and frequencies and acceleration levels identified. 

Early in the second year work, a successful 
deterministic rendezvous guidance algorithm based 
on optimal control theory was developed. The 
algorithm was considered sufficiently important 
that, with agreement of NASA, more emphasis was 
placed on the rendezvous investigation. To 
accommodate this work, the harpoon study was set 
aside. An effort was initiated on the rendezvous 
navigation problem wherein measurements are made 
and statistically processed onboard the spacecraft 
to provide the relative state information required 


for input to the guidance algorithm. Expenditures 
on the contract were low enough so that in June, 

1972, a no-cost extension of the work to September, 

1973, was possible. At this time the changes in 
objective were formalized and principal tasks were 
restated to Include the navigation work (and elim- 
inate the contact-rendezvous and harpoon investiga- 
tion) . 

In September, 1973, delays caused by installa- 
tion of a new computing machine at the University 
prevented generation of final data. The contract 
completion date was again extended, at no cost, 
to December 15, 1973, to allow time for this data 
generation and report preparation. 

The final report of our work is presented in 
two volumes: 

Part I, Guidance and Navigation Studies 

Part II. Tumbling Problem Studies 

Each of these volumes presents the technical 
details of the analyses conducted, the principal 
conclusions made, and listing of the computer 
programs employed, including descriptions of the 
operation of the programs. Technical abstracts of 
the work are included in each volume . In Part I , 
the body of the report reproduces a paper prepared 
for the AIAA 10th Electric Propulsion Conference 
entitled "Solar Electric Propulsion for Terminal 
Flight to Rendezvous with Comets and Asteroids." 
(AIAA Paper No. 73-1062). The title was changed 
for inclusion in thiB report and a few typograph- 
ical errors were corrected. 
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STUDY OF EFFECTS OF UNCERTAINTIES ON COMET AND 
ASTEROID ENCOUNTER AND CONTACT 
GUIDANCE REQUIREMENTS 


PART II. TUMBLING PROBLEM STUDIES 


Abstract 

The problem of determining the rotational motion 
of a tumbling celestial body of the asteroid type 
using space craft -acquired data is addressed. The 
rotational motion of the body is modeled by free- 
Eulerian motion of a triaxial, rigid body and its 
translational motion with respect to a nonrotating, 
observing spacecraft, which is not thrusting, is 
assumed to be uniform during the time observations 
are made. The mathematical details which form the 
basis for a digital simulation of the motion and 
observations are presented. Two algorithms for 
determining the motion from observations for the 
special case of uniform rotational motion are 
given. Results of studies of the performance of 
the algorithms using ideal data and "non-ideal" 
data are discussed. The "non-ideal" data is 
generated by inducing initial conditions which 
cause the body to "wobble" in general free- 
Eulerian motion. One algorithm performs satisfac- 
torily when the rotational motion of the body is 
truly uniform, but, in general, fails to determine 
realistic values of the constants which are deter- 
mined when the rotational motion is not uniform. 

The other algorithm does not perform well. 

I, Introduction 

Future missions to the asteroid belt or a comet 
using autonomous and/or semi -autonomous spacecraft 
of the CARD (Comet and Asteroid Rendezvous and 
Docking) type will foreseeably require the space- 
craft to land on, or "dock" with, the tumbling 
target body. The successful completion of such a 
task will depend to a large extent on the avail- 
able knowledge concerning the target’s rotational 
motion and its translational motion relative to the 
rendezvous spacecraft. Such knowledge, if not 
available a, priori , must be determined in some 
manner using observations of the motion of the tar- 
get body made by sensors onboard the spacecraft. 
Thus, methods for determining the rotational and 
relative translational motions of the target need 
to be developed. The development of such methods, 
or algorithms, is the subject of this report. 

A prerequisite of any motion 'determination al- 
gorithm, such as, for example, a satellite orbit 
determination algorithm, is the formulation of a 
suitable mathematical model for the motion which is 
to be determined. That is, the motion to be 
studied must be defined in physical terms and then 
described mathematically. The motion of interest 
here is that embodied in what we call "the tumbling 
problem." The definition of this problem is given 
the following section. Stated briefly, the motion 
consists of free-Eulerian motion of the target body 
about its center of mass and uniform translational 


motion of the body’s center of mass with respect to 
the center of mass of a nonrotating, observing 
spacecraft. Such motion can be described exactly 
In mathematical terms and the motion determination 
problem is reduced to the determination of the 
"constants of the motion" which are the arbitrary 
constants obtained in integrating the equations of 
motion of the system. 

Since the type of motion to be considered cannot 
be easily simulated physically, the mathematical 
model of the motion is useful not only in deriving 
the determination algorithms, but also in simulat- 
ing the motion digitally to check the performance 
of the algorithms. Such a simulation capability 
was developed and utilized in this study and is 
discussed in the third section of this report. 

The form of an algorithm for determining the 
constants of the motion depends, to a large degree, 
on the type of observational data available. In 
the fourth section of this document, two algorithms 
for determining these constants for the special 
case of uniform rotational motion {which will 
"probably" be the type of motion actually encoun- 
tered) are described. The first algorithm ,DUMRA 
(Determination of the Uniform Motion using Range 
and Angles), utilizes position data for three 
particular points on the surface of the target body 
at three instants of time. This data could be ac- 
quired by using the combination of (l) a scanning 
radar from which range and angles for many points 
on the body could be obtained, (2) a television 
camera, and (3) an earthbound observer. The 
observer would be used to identify the same three 
points on the body at the three instants of time 
and select the data for these points from among all 
that acquired (using, for example, a grided CRT 
display and a light pen) . 

The second algorithm utilizes simpler data. The 
data is range and range-rate data, also acquired 
via radar, for one specific point on the body at 
six instants of time, again with the aid of an 
observer to identify the point. 

The performance of the algorithms was tested, as 
indicated previously, by simulating the motion of 
the body and from this motion determining the 
observations required for the algorithms. The 
section following the description of the algorithms 
is devoted to a discussion of these performance 
studies. 

Conclusions drawn from the work conducted are 
presented at the close of the report and an 
Appendix containing information on the computer 
programs written and used during the study is also 
included. 
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II. Tumbling Problem Definition 

A prerequisite to the development of algorithms 
for determining the rotational state of a tumbling 
body, such as an asteroid, is the selection of 
suitable mathematical models for the body and its 
motion. Since the type of body of interest is an 
asteroid or a comet nucleus , it is reasonable to 
assume that the body is a rigid one. Also, using 
available data( 1 ) on the shape of asteroids, a 
triaxial, homogeneous, ellipsoidal, rigid body was 
selected as the physical model for the tumbling 
body. Of course, it is not suggested that any 
asteroid is actually perfectly homogeneous, nor 
that they all are ellipsoidal, but the shape ^ 
appears reasonable and the torque-free rotational 
dynamics of any rigid body can be taken as those of 
a body with a triaxial inertia ellipsoid. Our 
physical ellipsoid is the reciprocal of the cor-^ 
responding inertia ellipsoid and the lengths of its 
axes are denoted by a, b, and c where a >_ b > c. 

On the basis of observations of fluctuations of 
of light reflected by asteroids, astronomers ( 2 ) 
have concluded that they rotate uniformly, or at 
least almost uniformly, about their centers of mass. 
Flirt hermo re 5 the well-known fact that dissipation 
of rotational kinetic energy (internally, or 
through external means such as meteor impacts , with 
no significant change in rotational angular momen- 
tum will cause a triaxial body to ultimately rotate 
uniformly about its axis of maximum moment of 
inertia leads us to assume uniform, or almost uni- 
form, rotational motion of the tumbling body about 
its axis of maximum moment of inertia. To allow 
for nonuniform motion, we have adopted a free— 
Eulerian dynamical model for the tumbling motion; 
that is, the rotational motion is assumed to be 
torque-free during the time of observation. 

The relative motion of the spacecraft and the 
tumbling body centers of mass may be very nonuni- 
form during the use of high thrust chemical systems. 
However, we assume that during the time that obser- 
vations of the tumbling body are made the space- 
craft's high thrust system will not be operating. 

Of course, relative accelerations of the tumbling 
body and the spacecraft may still occur due to 
(l) the body's gravitational attraction, (2) the 
use of the spacecraft's low thrust system, and 
(3) the gravity-gradient effect of the sun's 
gravitational field. However, by assuming a sta- 
tion keeping or a steady fly-by mode of operation 
for the spacecraft., the relative acceleration 
caused by the low thrust system should be small 
enough to be neglected over a period of time of 
several minutes and the same should be true of 
the acceleration due to the body's gravitational 
attraction and also that due to the sun. Hence, it 
is assumed that the relative motion is uniform. 

The tumbling problem, as we have defined it, 
may thus be stated as follows: 

Using observations of the motion of a free, 
homogeneous, triaxial, rigid ellipsoid made 
from a reference point which is moving uni- 
formly relative to the ellipsoid's center 
of mass , determine the motion of the ellip- 
soid; i.e., determine the constants of the 
motion and the physical characteristics of 


the body needed to predict its relative 
position and its orientation at any time 
in the future. 

This is the problem we have addressed in this 
study. For the most part, however, it has been 
simplified by assuming uniform rotation as well as 
uniform relative translation. In the next section 
we consider the problem of simulating the observ- 
ations which a spacecraft might make. 


III. Simulation of Observations 

Although the algorithms for determining the 
tumbling body motion which are described in the 
next section are predicated on the assumptions of 
uniform translational and uniform rotational 
motions, simulation capability to handle nonuniform 
free-Eulerian motion of the triaxial body was de- 
veloped. Such capability may be useful in future 
studies and may also be used, as it was here, to 
study the effects of small variations in the rota- 
tional motion on the uniform motion algorithms' 
performance. Details of how the observations of 
points on the surface of the ellipsoidal model 
are generated using a FORTRAN subroutine, STATE 
(see Appendix for listing) are given in what 
follows . 

Physical Model 

The tumbling body model discussed in the pre- 
vious section is depicted in Figure 1. As stated, 
the model is a triaxial, rigid, homogeneous 
ellipsoid with axes of length a, b and c (a£b>c) 



Fig. 1. Tumbling Trolled C-eor.etry 
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so that the moment of inertia about the z-axis 
(see Figure l) is maximum. The moments of inertia 
about the axes x, y, and z are denoted by A, B, 
and C, respectively, and, as is obvious, C > B > A, 
in general O For a homogeneous, triaxial ellip- 
soid we have, 

A = (m/5) (b 2 + c 2 ) 

B = (m/5)(a 2 + c 2 ) (l) 

C = (m/5)(a 2 + b 2 ) , 


as shown in Figure 2 . The orientation of the 
Cx H y H z H system is defined by the angles % and 0 h- 
Finally, let the Cxyz system, shown in Figure 2, be 
a principal, body-fixed coordinate system for the 
ellipsoid. The orientation of the Cxyz system with 
respect to the CxjjyjjZjj is defined by the set of 
Euleri an angles {ij;, 6 , 4 1 ) ■ 

The rotational motion of the ellipsoid is 
governed by Euler's equations for a free triaxial 
rigid body, which may be written in the forms, 


where m is the mass of the ellipsoid. As we shall 
see, the mass of the ellipsoid does not affect its 
rotational motion as long as there are no external 
torques . 

Rotational Motion 

Let Sx^y^z^ be an orthogonal, dextral , coordi- 
nate system which has its origin, S, at the center 
of mass of the spacecraft and which is nonrotating 
Let CXYZ be a similar system which has its origin 
at C, the center of mass of the tumbling body and 

which is aligned with the 3x y z system. Also, 

S 3 S 

let Cx^y^z^ denote a similar system which has its 

z^-axis aligned with the constant rotational 
n 

angular momentum vector, H, of the tumbling body 


Y 



z 


Aa> x + (C-B) = 0 

Bui + (A-C) w w = 0 (2) 

y z x 

Cm + (B-A) oi w =0, 
z x y 

where “ x , uiy and ut 2 are the x-, y- and z - component s , 
respectively, of w, the angular velocity of the 
ellipsoid. 

Since Eqs . (2) are homogeneous, we can and shall 
divide each of them by a constant factor. In 
addition, we shall introduce a new independent 
variable , 

T = w a t » (3) 

o 

where w_ is the value of us z at t = 0« Letting 

z o 

O' = d{ ) /dx , we may then rewrite Eqs. (2) in the 
dimen s i onl es s forms , 


X' + (l-oO XX =0 
x 2 y z 

*V - (“r 1 ’ X A - 0 

*' Z ♦ (Oj-op Vy= 0 

where X /oj , X =oi /u , X =w /m , a = A/C 
x x z ’ y y z z z' z 1 
o o o 

and = B/C. 

Equations (4) have the first integrals, 

h 2 = 0 2 X 2 + cr 2 X 2 + X 2 = constant 
lx 2 y z 


a = (a, X 2 + a_X 2 + X 2 ) = constant, 
lx 2 y z 


where h = |h|/C 2 u> 2 and a is twice the rotational 
o 

kinetic energy of the body divided by C 2 ui 2 . The 

Z o 

integrals (5) and (6) may be used to obtain the 
following solution^) to Eqs. (4): 


X = p cn u 
x 

X y = q sn u (T) 


Fig, 2. Coordinate Systems. 


X = r dn u , 
z * 
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where 


p = 

{(0 - h^/tyi-o ^]} 55 

(8) 

q = 

{(a - h 2 )/[<J 2 (l-a 2 )]} !s 

(9) 

r = 

{(h 2 - cr 1 o)/(l-o 1 )]} !s 

(10) 

u = 

At - v, v — constant, 

(11) 

°2 ' 

- o 1 )(a-b 2 )/[(l-a 2 )(h 2 -d 1 c t ))} % 

(12) 


where A = 1 - (k 2 X^/q 2 ) sn 2 Xt. For use in 

Eqs. (l6) , p, q, r, X and k are to he obtained by 
evaluating Eqs. (8)-(l0)j (12) and (13) at t-0 
using Eqs. (5) and (6). 

The angles 6 and <j> may be expressed as functions 
of time by using the fact that H may be written as 

H = Aui x i + + Cui^k , (l8) 


or as 


and the functions cn u, sn u and dn u are Jacobian 
elliptic functions of modulus 

k = {[(o 2 -c 1 ) (a-h 2 ) ]/[(l-o 2 )(h 2 -o 1 a )]> 5 , (13) 

Note that, since Op and o 2 do not contain the mass, 
m, of the body, the rotational motion of the body 
is not affected by the value of m, but is affected 
by the distribution of the mass. 

To allow for arbitrary initial conditions on 
X x and X (initially X z = 1, of course), the 
elliptic functions are rewritten in the forms,( 4 ) 


cn u = cn (Xt-v) = cn Xt cn v + 
sn Xt sn v dn Xt dn v)/A 
sn u = sn (Xt-v) = (sn Xt cn v dn v - 
sn v cn Xt dn Xt)/A 


( 1*0 


H = - H sin 6 i + H sin 4> cos 4) j + 

H cos <(> cos 6 k, ( 19 ) 

where {£, j, £} is a set of unit vectors associated 
with the Cxyz system and H = |h| . From Eqs. (10) 
and ( 19 ) and previous definitions, it follows that 

sin 6 = - o^/h ( 20 ) 

and 

tan <f = CT 2 X y/ A z * ( 21 ) 

Expressing the angle 4* as a function of time is 
a more difficult task. Although it is well known( d ) 
that 4 / may be expressed as a linear function of 
time plus a periodic function of time composed ol 
Jacobi's Theta functions, this was not done. In- 
stead, 4 ) is calculated by integrating the follow- 
ing differential equation: 

4 >* = h [l +(a 2 -l) sin 2 <f> ] . ( 22 ) 


dn u = dn (Xt-u) = (dn Xt dn v + 

k 2 sn Xt cn Xt sn v cn v)/A , 


where 


A = 1 - k 2 sn 2 Xt sn 2 v. (l?) 

Since at t = 0 , cn u = cn v, sn u = - sn v and 
dn u = dn V, ve have from Eqs . ( 7 ) 

cn v = X x /p 
o 

sn \> = X /q (l6) 

y o 

dn v - l/r . 

Thus, the solutions for X x , \y and X z may be ex- 
pressed as 

A = (X cn At - A p/(qr) dn Xt sn Xt)/A 
x X o y o 

X = (X q/(pr) sn Xt + X cn Xt dn Xt)/A ( IT ) 

y x 0 y o 

X = (dn Xt - [X X k 2 /(pq) ] sn Xt cn At) /A 

2 y n 


Note that if the rotational motion is uniform, 

= 0 and 4> = Tp +hi . In the simulation program, a 
check is made to determine if the motion is uniform, 
and if so, the integration of Eq. (22) is by-passed 
and the orientation at the desired "time," t, is 
computed using 6 = 4 = 0 and 4> = 4* + hr . 

The transformation from a nonrotating set of 
unit vectors {I, J, K} directed along the positive 
Xs-, y s - and z s -axes, respectively, to the set of 
unit vectors 

(l j £) T = <i J k) t a t , ( 23 ) 


where a superscript 


A = 


10 0 
0 c$ s 4> 
0 -s$ c<j> 


c9 0 
0 1 
s 8 0 


T denotes the transpose and 


-s 6 


c 4 < s 4 ) 0 

0 


-s'|) c4 0 

c 6 


0 0 1 

_ 


_ 



c9 h 

0 

- S0 H 


a *H 

S *H 

0 

X 

0 

1 

0 


- S % 

o*„ 

0 


s8 h 

0 

C 0 H_ 


0 

0 

1 






_ 


(2k) 


In Eq. ( 2 U), c denotes the cosine and s denotes the 
sine of the angle prefixed by such letter. Since 
41 and 0 ^ are constant angles which specify the 
orientation of H, they constitute two of the con- 



stants of the motion. The other constants+are 

, X ,A , id and ^ Thus , we see that 
12x’y’z r o ’ 

O O 0 

eight constants are needed to completely specify 
the free-Eulerian rotational motion of a particular 
body . 

Translational Motion 

The model adopted for the relative motion of the 
centers of mass of the tody and the spacecraft is 
exceedingly simple. We assume that this motion is 
uni form so that the vector R from the point S to C 
is given by 

R = R o + Vt , (25) 

where 


The vectors rj smd pj may be expressed in the 
vector forms , 


ii = x j £ + y j J " + 

k 

(27) 

£j - 1 ♦ V ♦ ‘j 

K , 

( 28 ) 


so that in matrix form. 



5o = 


X I + Y_J + Z_K = constant vector 
o o o 


and 

V = XI + YJ + ZK = constant vector, 

are the relative position and velocity vectors at 
t = 0. The constant components of R and V are 
six constants which determine the relative trans- 
lational motion. 

Motion of a Point on the Body 

The determination algorithms discussed in the 
next section are predicated on the use of observa- 
tions of the movements of specific points on the 
surface of the tumbling body. If we let rj be a 
vector from C to a point , F j , on the surface of the 
body and pj be a vector from S to Fj (see Figure l) 
then we have, very simply. 




= R + r. 


(26) 


(h 


where 



(30) 



-A A 0 
y x 


Equations (29) and (3l) are used in subroutine, 
STATE, to generate the relative position, p j , 
and relative velocity, p j , of three points, Fj , 
j=l,2,3, on the surface~of the ellipsoid. The 
range, P j , and range-rate, P j , for each point are 
also determined from 


P 


j 




+ n 


+ 5 2) ' 


(32) 



and 

(33) 

respectively . 

As stated above, the points P^ are "on the 
surface" of the ellipsoid. This constraint is 
satisfied by allowing one to input only the x- 
and y- components of rj and determining the 2 - 
component from 

Zj - c (l - x 2 /a 2 - y^/b 2 )* 2 ^ < 3h ) 

where a, b and c are to be provided and x^ and y^j 
are constrained by the relations xj £ a, yj _< b, 

1 - xj 2 /a 2 - yj 2 /b 2 >_0. Note that we have assumed 
that the points are all "on the bottom" (i.e. , 
z >_ 0) of the ellipsoid. 

Observability of a Point on the Ellipsoid 


Fig. 3. Geometry for Visibility Determination. 


t Actually, and are physical parameters. 


Clearly, a point Pj , on the ellipsoid will not, 
in general, be visible from a point S at all 
times. To determine whether or not Pj is visible, 
it is sufficient to find the two (in general) 
intersections of the line of sight, Lj , which is 
colinear with p j , with the ellipsoid (see Figure 3) 
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and determine if the point, Pj, is the near point 
of intersection . Such a determination is made in 
subroutines STATE and VISIB in the following manner. 

The components of the vector, p j , in the {i,j,k} 
basis are determined using 


Using the coordinates (x*, y* , z*) the distance, 

J ti J 

p*, may be computed and if p* > P., the point, P, 
J j J J 

is visible, but if pj > P*, the point, Pj , is on 

the "backside" of the ellipsoid and cannot be 
observed. 


* ' 
A 


1 

e j 



A 



n J 

11 

11 > 

n J [ 

A 



C J 



1 4 


t * 


(35) 


and the direction cosines, 2.. , ra. and n. , of L. 

J J J J 

are determined using the equations. 


= «/ /P J 


m Aj " n j /p ] 


(36) 


"Aj * b /p J 

« j, # 

The coordinates x^ , y^ and Zj of the "other" 

point of intersection of Lj with the ellipsoid must 
satisty the equations, 


* - x ) / £ = (y* - y,)/m, = (z* - zj/n 


J J 


A J j 


a" a 


j y f j 


and 


(37) 


1 = (x*/a) 2 + (y*/b) 2 + (z*/c) 2 . 

<} 0 v 


. t 


The solution of Eqs. (37) and (38) is 


(38) 


X* = 2 o 2 /c 1 - Xj 

= <m A. /e A J !<X * - x j> + y J 

z j = (n A - V + V 


where 


(39) 


-1 = £ A 2 /a 2 + m 2 /b 2 + n 2 /c 2 

1 j j J 


(1*0) 

2 


c 2 = (o l * V /a *> X J - m A , W” 2 ' "a/A-V" 


J A 


A A 


Ue note that if Z = 0 the second and third of 

J 

Eqs. (39) appear to be singular. This singularity 
is, however, only "apparent," since 

lim (x*-x.) ~ 2[y j {m A /b2) + z j (n A /c2)] 

£ -> o f = J J 

a 3 [m Aj /b +n Aj /c2] (4!) 


In summary, the motions of three points on the 
ellipsoid are simulated digitally in subroutine, 
STATE, which utilizes Eqs. (29) and (30) (and all 
equations needed to evaluate those equations). 
Values of a, b and c, as well as R q , V, 6^, 

w = U v i + \ r A + k) , x and y , J=l,2,3, 

"° Z o o y o J J 

must be supplied to STATE by the user. As the 
ellipsoid rotates and translates, the visibility 
of the chosen points is checked by using Eqs. (39), 
(32) and (29). If the point, Pj , is not visible, 
a message to that effect is printed, but, since 
subroutine, STATE, was used only to check out the 
motion determination algorithms, the data request- 
ed by the algorithms is still supplied. In future 
simulation work, however, the visibility provision 
should be used to make the simulation more real- 
istic by "cutting off" data from an unobservable 
point. 


V. Algorithms for the Case of Uniform Motion 

Because asteroids and comet ary nuclei may 
logically be expected to be rotating uniformly and 
because the exact solution to the free-Eulerian 
motion problem is much more complex if the motion 
is not uniform, algorithms for determining the 
eleven constants, , 9^, w = |m| , X, Y, Z, X q , 

Y , Z x x and y^ , were developed. For the uniform 

motion case,^ and determine the orientation of 
the target’s angular velocity vector w, which is, 
of course, colinear with its rotational angular 
momentum vector, H, and m is the magnitude of the 
singular velocity. The constants X, Y, Z, X 0 , Y q 
and Z Q determine the translational motion of a 
point, C, on the spin axis of the target and and 
y-|_ locate a point, Pq , with respect to C. The z- 
coordinate of is not determined, since the loca- 
tion of the center of mass of the target on the z- 
axis cannot be found if it does not "wobble." 

That is, if the motion is uniform, all points on 
the axis of rotation move in the same manner as 
the center of mass and hence are not distinguish- 
able from it. This fact is taken into account by 
setting = 0. Furthermore, since the orientation 
of the Cxyz system in the asteroid is immaterial if 
the motion is uniform, ij/ can be set equal to zero 

i 0 

also . 

Two distinct algorithms were developed for the 
uniform motion case. The first is called DUMRA 
(Determination of the Uniform Motion using Range 
and Angles). This algorithm is presently set up 
to accept position vectors , which it is assumed 
may be obtained from range-and-angle data for three 
specific points on the target's surface. The 
second algorithm is called DUMRRR (Determination 
of the Uniform Motion using Range and Range-Rate) . 
This algorithm utilizes range and range-rate data 


A quadratic equation in x* may be obtained. One root of this equation is obviously x so that, the 
first of Eqs. (39) may be J obtained by factoring the quadratic. J 


6 



for one point on the target's surface along with 
Estimates of the eleven constants to he determined 
&nd iteratively produces more exact values for the 
constants. More detailed descriptions of these 
algorithms are given in the remainder of this sec- 
tion. 

Algorithm DUMRA 


and t£ has been arbitrarily set equal to zero. 
Equation (1*8) may be solved for and y lt viz., 


4 £i • 

(5*0 


M . 

[2(cmAt - I)]’ 1 0 Q 

M 

0 [2(cu>At - l)] _1 0 


The solution procedure utilizing position data 
is as follows: 

1. Obtain position vectors, p j , j=l,2,3, for 
three distinct, non-colinear points, P j , j=l,2,3, 
on the surface of the target at three times, t j , 
j=l,2,3, such that t 2 - t-j_ = tg - t 2 - At. 

2. Compute approximate velocity vectors for 
each point at time, t 2 , using 

Pj(t 2 ) = [pj(t 3 ) - p J (t 1 )]/2At, j=l ,2 ,3. (U2) 

3. Compute the angular velocity of the body, 
using! b ) 

[p_(t_) - p At-)] x [p 9 (0 - p _ ( t 0 ) ] 
m =2_J 2 - 1 — 2 . , , (U3) 

£p 3 (t 2 ) - (p 2 (t a )] - (p 2 (t 2 ) - p x (t 2 )] 


so that in matrix form, 

m = [u 1 u> 2 vi^]' 

4. Compute 4>tr and 6™ using 

n n 


= tan -1 (m 2 /tA> 1 ) 


and 


where 


9 H = tan 1 ( /1-cos © H / C os0„) » 


cos 0 H = u 3 /|ta| . 

5. Compute x^ and y^. First, compute 

APl = - 2 P-^g) + £1^1^ 

in the matrix form - 


(WO 

(1*5) 

( 46 ) 

( 1 * 7 ) 

( 1 * 8 ) 


6. Compute R(t 2 ) = R . This may be done using 
the equation ~ 


5o = £l (t 2 } - ^3 ^2 <1 y l 


T. Compute V q from the equation, 


!o = [ £ ( V ]/At 


(55) 


(56) 


The requisite constants are given by Eqs. 

(1*3) » (45), (46), (54), (55) and (56). Note that 
to extract the coordinates and yi, it was 
necessary to allow for the curvature of the motion 
of due to rotation. 

Algorithm DUMRRR 

The solution procedure based on the use of range 
and range-rate data is an iterative one. The main 
steps are as follows: 

1, Guess values for the eleven unknowns, 

0 R , w, X q , Y qJ Z , X, Y, Z, Xl and y r 

2. Set ^ = 0 at t = 0 and compute, for six 

times t >0, 

J 

p, = (R T R + 2 R T A r, + r? rj„ (57) 


- ail ' 11 -1 e 


and 


0 rn rprp** rp rp 

( P P-i ) = (V; R + R 1 A 1 u r + V*" A 1 r ) 

lie -0 - — = = — 1 -o = -1 e 


(58) 


where 


.T .T 
A 2 

[*' 

T (t 3 ) 

-21+ A T (t 2 )] 

y i 

. (49) 

where the subscript e 

denotes an 

e st imat ed value , 







0 



"0 -rn 0“ 






cot scat 

0 




m - 

Ci) 0 0 


(59) 

4i (t) 

= 

-suit cot 

0 



(50) 


0 0 0_ 






0 0 

1 




T T T T 

and A = A 3 A 2 A£ . 




I = 

3 

x 3 identity matrix 

j 


(51) 

3. Use the observations p = 

|p^| and 

p f p i /p i 



C0 H 

0 - S 0 H - 





at the same times, t., 

and the estimated 

values to 

=2 


0 

1 0 




(52) 

obtain J 

( 

\ 




- S0 H 

0 c0jj 





Az = 

/>-; 

le 

( 60 ) 



c*„ 

"*H 0 






( p i°i - < 

p l p l>e) 


&3 = 


-S *H 

C *H 0 

y 



(53) 

at each time, t j , j=l,: 

2,... ,6. 





0 

0 1 
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!+, Evaluate the partial derivative matrix. 



where X = (t^ 0^ id Y q X Y Z y^) , at each 
time, t j , j=l,2,...,6. 

5. Form the augmented matrices. 


and compute Ax from the equation, 

Ax = D g T Az , ( 6 7) 

at successive instants of time, correcting the 
estimate, x e , each time. This method of solution was 
attempt ed"initially when there were still errors in 
the program. Since it appeared that the matrix 
H TH was ill-conditioned, the pseudo-inverse method 
was = abandoned in favor of the more cumbersome, hut 
more trustworthy, "batch-processor" method outlined 
here . 



(62) 


VI . Performance of the Algorithms 


The algorithms, DUMRA and DUMRRE, were tested 
using the simulation capability provided by sub- 
routine, STATE, (and subsidiary subroutines). They 
were not exercised as much as was desired because 
of time constraints, hut whether they operate satis- 
factorily in the uniform motion cases considered is 
evident from the results discussed in^this section. 
Nonuniform motion due to target body "wobble" is 
another matter and is also discussed infra. 


and 

C= [^ T (t 1 );'^\t 2 K . .;'Az rr (t 5 ); (p x - Pl ) tBt ] T 

(63) 


To test the performance of the algorithms , 
values of the physical parameters a, b and c 
thought to be .typical of an asteroid were chosen. 
These values are given in Table 1 along with 
values of the nine constants of the motion and the 


6. Invert B and compute 

Ax = E,- 1 C C f 


TABLE 1 

PHYSICAL PARAMETERS AND CONSTANTS OF THE MOTION 


where C f is a convergence factor (C f £ l) which may 
be varied to improve convergence. 

7. Find a new estimate of x using 

X i+1 = X 1 + 4* . (65) 

-e -e — 

8. Check the norm, | | C_ j j > to see if it is suf- 

ficiently small. If it is, go to 10. 

9. To decrease | [ C_[ | , compute a new partial 
derivative matrix B and a new matrix C using the 
current estimates and repeat steps 6 through 8. 

10. The elements of x e 1+1 are the desired 
values of the constants. 

The matrix § in Eq. (62) is such that some of it 
Its terms are much larger than others. Therefore, 
to increase numerical accuracy in obtaining the 
inverse of §, dimensionless variables are used in 
the computer subroutine DERIV in which the matrices 
H(tj) are computed. 

It should be noted that the method of solution 
used in this algorithm is of the "batch processor" 
type. Another method for solving such non-linear 
observation equations as Eqs. (57) .and (58) is to 
use the pseudo-inverse, 

a = (a T gr 1 , ( 66 ) 


a = 1.75 x 10 4 m w = 3.31 x 10 -4 rad/sec 

b = 8.0 x 10 3 m 0 H = 20° 

c = 3.5 x 10 3 m i|> H = 10° 

X = - 1.0 x 10 6 m u = 1.0 m/sec 

o 

Y = 0.0 v - .5 m/sec 

o 

Z = 0.0 . . w = .5 m/sec 

o 

x = - 1.75 x 10 4 m y 1 = 0.0 

x 0 = - 1.0 x 10 4 m y 2 - 6.5652 x 10 4 m 

x = - 1.0 X 10 4 m y = 0.0 

3 ^ 


X- and y-coordinates of three points on the body . 
The z-coordinates of these points are determined 
in STATE as indicated previously. The values of 
the constants which are listed constitute only one 
of several .sets of values used during the check- 
out process and hence are "typical." 

The first step in the check-out of each of the 
two algorithms was to test its. accuracy under 
ideal conditions represented by perfectly uniform 
motion. The performance of each algorithm when 
the simulated motion was not uniform was then 
tested. 
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For perfectly uniform motion, algorithm, DUMRA, 
reproduced the input in most cases to four signify 
icant fi gures and always attained percent errors 
less than 1 percent. Theoretically, the only error 
which is present (apart from basic computational 
errors, that is, is that due to the approximate 
calculation of the velocity of each of the points 
from the position data. This error does not appear 
sufficiently large to be a problem if the increment 
in time. At, is small compared to the rotational 
period of the target. In the cases examined here, 
u> = 3-31 x 10“ 4 rad/sec, so that the period of 
rotation is 1.9 x 10 4 sec = 5.28 hrs . For these 
numerical values. At = 90 sec., gave satisfactory 
results . 

When nonuniformity of the rotational motion was 
introduced by giving w Xo and/or u non-zero values , 
the effect on algorithm DUMRA was' anomalous when the 
maximum value of the angle, O = cos^Ccos 0 cos 4 *), 
between the rotational angular momentum vector and 
the angular velocity vector was greater than a few 
degrees. The algorithm extracted fairly accurate 
(within 10 percent) values for the components of R 
and V c at times, but failed as often as it was suc- 
cessful. The reason for this behavior has not been 
determined at this time, but appears to be connected 
with the use of the position vectors to calculate 
velocity vectors. 

The performance of algorithm, DUMRRR, was not 
satisfactory. In the runs made, the use of the 
algorithm usually resulted in the reduction of 
square root of the sum of the squares of the ele- 
ments of the matrix, Q, but this reduction did not 
lead to consistent improvements in all the constants 
which must be determined. The errors in R 0 were, 
for example, reduced by applying the algorithm, but 
those in V c and r, in general were not. In a few 
cases, application of the algorithm resulted in 
progressively smaller values of j |c| | for several 
iterations, but continued applications resulted in 
divergence. 

Attempts were made to improve the convergence 
properties of the algorithm by varying the conver- 
gence factor, Cf. An initial value of 0.02 was 
given to Of and this value increased to 0.1 and 1.0 
when convergence criterion was reduced to 10“ 4 and 
10 -6 , respectively. These efforts were not suc- 
cessful . 

Since the algorithm, DUMRRR, did not perform 
satisfactorily when confronted with uniform motion, 
tests with nonuniform rotational motion were not 
conducted. 


Conclusions and Recommendations 

In summary, the capability for simulating the 
motion of an asteroidal-type body and observations 
of points on such a body has been developed. The 
rotational motion simulated may be general free- 
Eulerian or uniform, while the translational motion 
is uniform. This simulation capability has been 
used to test two algorithms for determining the 
constants of the motion from simulated observations 
when the rotational motion is uniform. 

The algorithm, DUMRA , performed the job for 
which it was designed satisfactorily. However, 
when required to work in the presence of "noise” 
caused by nonuniform rotational motion, the algo- 


rithm produced erroneous results. Failure of the 
algorithm for large nutation angles (angle e) was 
expected, but even small angles unexpectedly caused 
a great deal of trouble. Because, small nutation 
angles may actually be encountered, the algorithm 
should be modified to allow for such an eventuality. 
This could be done by using the angular velocity 
determination portion of the algorithm with a small 
At to obtain five values of 10 . An approximate 
solution which allows for small nutation angles (see 
Ref, 5, Appendix D) could then be used to obtain 
w x , uiy , and c^. Minor changes in DUMRA 

would l Hen a?low the center of mass location and 
motion to be determined. 

The work on algorithm, DUMRRR, based on the use 
of range and range-rate data, although not success- 
ful, should not be considered to have been fruit- 
less. The failure of this algorithm reinforced an 
intuitive belief that more information than just 
range and range-rate data are needed to determine 
the motion. In particular, it appears that angular 
data is necessary to adequately define the plane of 
rotation of the vector and hence the direction 
of to. It was hoped that the requirement for such 
data could be eliminated by having very good esti- 
mates of the angles and 0j[, but this does not 
seem to be the case. 

An important aspect of this study is that an 
observer was assumed to be available for the 
purpose of identifying points and selecting data. 
Totally autonomous operation of course precludes 
the use of such an observer . Thus , some type of 
automatic identification technique should be de- 
veloped if a totally autonomous determination of 
the target body’s motion is a requirement of 
actual missions. 
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APPENDIX - COMPUTER PROGRAMS 



The DUMRA and DUMRRR programs are written in 
FORTRAN IV and executed on the IBM 370/155- The 
programs are research tools, not developed produc- 
tion routines. The steps in the simulation, and the 
names of the subroutines used, are briefly as 
follows. 


DUMRA 


Values of time, T, and the time increment, DT, 
are input to subroutine, MAIN, which calls sub- 
routine, DUMRA, which calls subroutine, STATE, at 
times, T1=T , T2=T + DT and T3=T + 2*DT. STATE 
receives input from subroutine, DATA, the first 
time it is called and produces the position vector 
sets, R1 (I,J), R2 (I,J), R3 (l,J), I = 1,3, J = 1, 
3, for times, Tl, T2 and T3, respectively. Here, I 
is the coordinate index; i.e., I = 1, 2 or 3 im- 
plies the x-, y- or s-coordinate , respectively, of 
a point, P j . The integer J is the number of the 
point. The position vector sets are used to de- 
termine the constants of the motion. In both STATE 
and DUMRA matrix algebra subroutines are used to 
perform the necessary calculations. Subroutine 
STATE, also calls DJELF and DCEL1 during the simu- 
lation process. 

DUMRRR 

The MAIN subroutine of this program executes the 
DUMRRR algorithm described in the body of this 
report. Subroutine, MAIN, accepts a time, T, 
a time increment, DT, and estimates of the con- 
stants of the motion. It calls subroutine, RHOES, 
which in turn calls STATE and computes a difference 
matrix DZ (l,l), I = 1,2, formed from the differ- 
ences in the actual range, RHO , and range times 
range-rate, and the estimated values of range, ROE, 
and range times range-rate, R0R0DE, respectively, 
for the times at which it is called. Subroutine, 
DEHIV, is also called from MAIN for each of the 
six times necessary and computes the matrix H(I,J), 

I = 1,2,..., 11, J = 1,2 of partial derivatives of 
the observations with respect to the constants being 
determined. Matrix algebra subroutines are called 
from MAIN, RHOES, STATE and DERIV as needed. The 
observations are stacked in an 11 x 1 vector, 

C(l,l), and the partial derivatives are stacked in 
the 11 x 11 matrix,B(l ,J) . Increments in the 
constants of the motion are calculated as an 11 x 1 
vector, DX(l,l), and new estimates obtained. The 
algorithm is repeatedly executed until the norm 
of the differences in the actual and estimated 
observations is "sufficiently" small. 



t 


Subroutines are used in both DUMRA and DUMRRR. 


STATE Simulates the free-Eulerian Rotational 
motion and the uniform translational 
variation, produces observational data 
and checks the visibility of the ob- 
served points. 

DUMRA Determines the eleven constants, PSIH, 

THETAH, WW R , V , x.. and y, . 

— o — o ± l 

MATMUL Matrix multiplication subroutine. 

TILDE Matrix algebra subroutine which forms 

a skew-symmetric, 3 x 3 matrix from a 
3x1 matrix (vector) for computation 
of a cross-product. 

D0TPR0 Matrix algebra subroutine which forms 

the dot product of two vectors . 

DATA Reads in data needed in STATE. 

VISIB Checks data from STATE and prints out a 
message if any of the three points being 
observed is (are) not visible. 

DJELF Computes the Jacobian elliptic functions 

snu, cnu and dnu. 

CELI Computes the complete elliptic integral 

of the first kind of modulus k. 

Entries are indicated in the flowchart 
in Figure A-l. 

DUMRRR 

MAIN Controls execution and executes the al- 
gorithm, DUMRRR. 

STATE Same as for DUMRA except for minor 
modi fi cations. 

RHOES Computes the difference Ajs in the actual 
and estimated observations. 

DOTPR Same as for DUMRA. 

TRPOSE Transposes an input matrix. 

TILDE See DUMRA above. 

MATMUL See DUMRA above, 

DATA Reads in actual data. 

A37I^ Inverts matrix B. 

DJELF See DUMRA above . 

DECL1 See DUMRA above. 

Entries are indicated in the flowchart . 
in Figure A-2. 
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•Called from DUMRA and STATE. 


Fig. A-l. Simplified Flowchart for DUMRA 


Controls execution and contains 


DUMRRR algorithm. 
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Variable Names and Definitions 

The most important variables in regard to INPUT 
and OUTPUT are defined in the sections which 
follow,- Other variable names are defined through- 
out this appendix. 


Input Data 

The data which must be provided as input to 
programs DUMRA and DUMRRR is described below. All 
field specifications are F16.0 and the card 
column in which a particular piece of data must be 
located are indicated by numbers in parentheses 
following the mathematical symbol for the piece of 
data. 

DUMRA 

Card 

No. Information on Card 

1 Initial time and time increment: t c (ccl-l6), 

At ( cc 17-32), 

2 Constants of the Motion: (cc l-l6), 

to (cc 17-32) to (cc 33-46),° (ec 49-64), 
y z 0 

o o 

0 H (cc 65-80). 

3 Constants of the Motion : <1 >jj( cc l-l6) 

X(cc *+9-64), X Q (cc 65-80). 

4 • Constants of the Motion: Y (cc l-l6) , 

Z ( cc 17-32). 
o 

Coordinates of Points: x^(cc 49-64 ), 

y^(cc 65-80), x 2 (cc 65-80) 

5 Coordinates of Points: y 2 (cc l-l6), 

x^(cc 17-32), y^(cc 33-48). 


Card 

No . Information on Card 

6 Actual Values of: ^(cc l-l6) , x(cc 17-32) , 

y( cc 33-48), Z(cc 49-64), X q (cc 65-80) 

7 Actual Values of: Y q (cc l-l6) , Z Q (cc 17-32), 

x 1 (cc 33-48). y x (cc 49-64), x 2 (cc 65-80) 

8 Actual Values of: y 2 (cc l-l6) , x^(cc 17-32), 

y 3 (cc 33-48) 

9 Physical Parameters: a(cc l-l6) , b(ce 17-32), 

c (cc 33-48)) 

Output Descriptions 

The outputs of the programs, DUMRA and DUMRRR, 
are self-explanatory if one refers to following 
definitions of quantities which appear on the 
printouts. 

DUMRA 

WXO, WYO, WZO = initial values of the components of 
angular velocity in the Cxyz system, 

PSIO = initial value of , 

THTAH = 0 . 

PSIH = * . 

XD, YD, ZD = X, Y, Z, respectively. 

X, Y, Z = X q , Y , Z q , respectively. 

SX, SY, SZ = x-, y- and z-coordinates , respectively, 
of the point indicated. 

PHI = 4> - 
THETA 5 6 . 


6 Physical Parameters: a(cc-l-l6), b(cc 17-32), 

c(cc 33-48). (Semi-axes of ellipsoid) 

DUMRRR 

Card 

No. Information on Card 

1 Initial Time and Time Increment: TC(BC l-l6) , 

DT(cc 17-32) . 

2 Estimated Values of: tf'jjCcc l-l6) , 

e u (cc 17-33), to ( cc 33-48), X (cc 49-64), 
n o 

Y (cc 65-80) , 
o 

3 Estimated Values of: Z (cc l-l6) , x(cc 17-32) , 

o 

Y(ec 33-48), Z(cc 49-64), x^(cc 65-80) 4 


WW = |uj 

W = angular velocity velocity, w. The x y s ~ and 
z s -components of w are listed under W. 

V = V - The components, X, Y and Z, are listed 
under V. 

R0 = R(t 2 ), the vector R at time, TC + DT. The 
components of R are listed under R0. 

SRI, SR2 , SR3 = the components of r^» r 2 , and r 0 , 

respectively, in a coordinate system 
rotated with respect to the Cxyz system 
so that 1(1 = 0 at time, TC + DT, and 
translated so that = 0 are printed 
underneath these symbols. 


4 

5 


Estimated Value of: y^Ccc 1-16) 


Actual Values of: u (cc 1-16) , w (cc 17-32), errors 

x o y o 1 

u (cc 33-48), (cc 49-64), e u (cc 65-60) the components 

Z OH 

O 


actual-estimated 

actual , in the values for 

of R^t^), r^, r Q , r^ and V^ computed 
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using simulated values as true values, are printed 
also under the headings DRO, DSR 1, DSR 2, DSR 3, 
and DVO, respectively. Similar errors in the 
components of fi , r 2 and r^ are also given under 
the headings DRDOT 11, DRDOT 2 and DRDOT 3, re- 
spectively. It should be noted that when 0.999 D03 
is printed this does not indicate a 99,900 percent 
error, but is an indication that the actual value 
and the calculated values are both zero. 

DUMRRR 


The output of DUMRRR includes the values of the 
constants of the motion and the other data needed 
for the simulation of the motion in the same out- 
put format as that used in DUMRA. Additional out- 
put includes the estimated values of the constants 
of the motion and the coordinates of P^, the 
times, Tj at which data is taken, the vector C 
formed from the actual and estimated observations, 
the estimated values of the constants PSIH, THETAH, 

WW (WE on printout ),R , V and x, and y, and the 
-o -o 1 1 

matrix, DX, of changes in the constants. 
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non 


DUMRA Program Listing 


IMPLICIT REAL*8(A-H,0-Z> 
NN=l 

READ( 5, 1 )TC,DT 
1 FORMAT I 2F 16. 0 ) 

DO 100 1=1,100 
NC = 1 

CALL DUMRA (TC,DT,NC»NN) 
TC=TC+DT 
100 CONTINUE 
STOP 
ENO 


SUBROUTINE DUMRA ( T 1 , DT , NC , NN > 

IMPLICIT REAL*8(A-H,0-Z I 

DIMENSION VS (3, 1 ) , ROS I 3 , l ) , SRRS ( 3 , 3 > , ROOTS ( 3 , 3 > , OSR ( 3 ♦ 3 

• I, DV (3,1), DR 0(3,1), DRDOT (3*3) 

DIMENSION Rl(3,3) ,SR1(3,1), SR 2(3,1) , SR 3(3,1), ROOT <3,31 
DIMENSION WI (3) , SR I ( 3, 3 I ,RH0 1 3 1 ,RH00(3I ,R2(3,3> 
DIMENSION R3 ( 3, 3 ) « DR L ( 3 , 3 ) , DR2 ( 3 , 3 ) » RDDOT { 3 , 3 ) , W ( 3, 1 1 
DIMENSION DELRK3, 1) , H( 3, 3 ) »0ELR01 ( 3 , 1 ) , DELR02 (3,1) 
DIMENSION ATAP( 3*3) ,ATH(3,3J , DELT ( 3 , 3 I , VCV ( 3 , 1 1 ,V(3,l) 
DIMENSION WT ( 3 » 3 I , RD ( 3 , 1 ) ,V2(3,1),R0I3,1)*SR(3,3) 
DIMENSION SR 11(3,1), SR 12(3,1), SR 13(3,1) 

DIMENSION ATAPT(3,3) ,A3( 3,3) ,A3T(3,3) ,SR1T3(3,1) 
DIMENSION A3T1( 3, 3) , A3T3( 3 ,3) ,DATI (3,3) ,DDR1(3,1) 
DIMENSION E(3,3),AT3(3,3),AT1(3,3),XI(3) 

DIMENSION DRU(3.1),02R1(3»1) 

DIMENSION SRI L( 3, l ) 

DIMENSION SSR<3,3) ,VS1(3,1), R0S1 (3,1) 

DIMENSION AB(3,1) ,BC(3,3) 

DO 210 1=1,3 
00 209 J=l , 3 
R 1 ( I , J )=O.ODO 
R2( I , J ) =O.ODO 

209 R3 ( I , J ) =0 • ODO 
V2( 1,1 )-0«0D0 

210 V( 1,1) =0*000 

OBTAIN OATA 

CALL STATE (Rl, SR, SRI ,RDOT , WI , DPS I , T 1 , T , NC , I V I S , RHO, RHOD 
, ,TTT,PPP,NN, VS1 ,ROSl,SSR) 

CALL VISIBIIVIS) 

PRINT 997, T 
T2=T1+DT 
NC = 3 

CALL STATE (R2, SR, SRI , ROOT , W I , DPS I , T2 , T , NC , I V I S , RHO, RHOD 
« ,TTT , PPP,NN, VS1,R0SI,SSR) 

DO 116 1=1,3 
VS( I , 1 )=VS1 ( 1,1) 

ROS ( 1,1)=R0S 1(1,1) 

DO 116 J= 1 # 3 
SRRS( 1 , J ) = SSR ( I ,J) 

116 ROOTS (t,J)=RDOT( I ,J) 

PRINT 200 , DT , T , T2 

200 FORMAT (/, 5X» * DT = * , F12.6, 10X, * T = • , F 1 2 . 6 , 10X , • T2 = •, 

• F 12 » 6 « / I 

CALL VISIBIIVIS) 

PRINT 310, DPSI 
T 3=T 1 *2. 0*DT 
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CALL STATE I R3» SR, SRI , RDDT f W I * DPS I , T3, T,NC, I VI S * RHO, RHOD 
. ,TTT,PPP.NN,\/Sl,ROSl,SSR) 

CALL VISIBIIVISI 
PRINT 947, T 

947 FORMAT(5X,'T='»D16.8*/) 

3LO FORMAT! 10X, ' PSI = ',020*8,//) 

C 

C FIND RDOT AT TIME T2 

C 

DD 101 1=1,3 
DO 101 *1=1,3 
DRi ( I , J ) =0.0 
0R2( I , J)=0.0 
101 ROOT I I » J ) =0-0 
00 10 J= 1 * 3 
DO 10 1=1,3 

DR1U , J) = CR21I, J>-RKI,J» >/0T 
DR2( !,JI = !R3U, J)-R2(I, J> )/DT 

10 ROOT! I , J > = !DR1! 1 1 JM-DR2! l, J) )/2.0D0 
C 

C FIND W 

C 

IF(DT.GT. 100*0) GO TO 543 

DO 11 1=1,3 

OELRK 1,1) =0* ODO 

DELRD1 ( I,1)=0.0D0 

DELRD2 I I , 1 )=0.000 

OELRK I ,1)=R2( I ,2)-R2(I ,1) 

D6LRDK l,l)=ROOT( I ,2)-RD0T (1,1) 

DGLR02 ! I,i)=ROOTt |,3)-RD0T( 1,2) 

O0TP = 0 *000 
VCVI I,1)=0.000 

11 XUI )=DELRD2 (1,1) 

CALL TILDE ! X I »DELT I 

CALL MATMUL ( DELT,0ELRD1* VCV* 3*3* 1 ) 

CALL DOT PRO! DELRD2»DELRl ,DOTP) 

DO 15 1=1,3 

15 Ml I V U=VCVU , 1 ) /DOTP 

WS=W(1, t )*W! I,L)+W(2, l)*W(2»l>+W(3« l)*WC3,l> 

C 

c FORM APSI ANO ATH 

C 

543 WW=DSQRT(W$! 

PRINT 544, MM 

544 FORMAT!/, 5X, 'MW = , ,D20.8,/J 
DSI=0$lN( HW*DT ) 

DCO=OCOS ! WW*DT ) 

CTH=WI3,1)/WW 
STH=OSQRT! 1.000-CTH4CTH) 

TH=DATAN( STH/CTH) 

DTH= TH# 180. 0/3. 14159 
IF!W<2,U ,LT.i.06-6)GO TO 27 
P$1H=DATAN2! Ml 2,1 ) ,M! 1 , 1 ) ) 

GO TO 28 

27 P S I H=0 • 0 

28 PSIHD=PSIH*180* 0/3. 141592 
CP=DCOS!PSIH) 

SP = DS I N ( PS I H ) 

DO 22 1=1,3 
00 22 J= 1 , 3 
ATH! I , J ) =0 . ODO 
22 HI I , J)=O.ODO 
ATH! 2,2)=1.0D0 
ATHll,l)=CTH 
ATH! 1,3) =- STH 
ATH! 3 , 1 ) =STH 
ATH(3,3)=CTH 
H 1 3 , 3 ) = 1 .000 
HI l, 1 )=CP 
H( 1 , 2 ) =SP 
H!2,n=“SP 
H(2,2 >=CP 

CALL MATMUL! ATH , H , AT AP , 3 . 3 , 3 ) 

DO 120 1=1,3 
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ono ooo o o o ooo ooo ooo nno on 


00 120 J = l,3 
A 3 ( I t J )=0. 0 

120 ATAPTH* J>=ATAP< J.I ) 

FORM A3T1 ANO A3T3 

00 30 1=1*3 
DO 30 J = l* 3 
E ( I , J ) =0*000 
A3T1 ( I * J ) =0*000 

30 A3T3( I* JJ=0.000 
A3T1 (1, 1)=DC0 
A3T3 ( I t l ) =0C0 
A3T1 ( L » 2 ) =0S I 
A3T3< 1,2)=-DSI 
A3T 1(2*1 )=-D$I 
A3T3< 2 * 1 ) =D$ I 
A3T1 ( 2 1 2 l = DC0 
A3T3( 2,2 }=0CQ 
A3T 1 ( 3 » 3 ) = 1. 000 
A3T3(3, 31=1.000 
00 31 I = i » 3 
DDRK It 11=0.000 
02R H 1,11-0* 000 

31 6U,n = 1.000 

FORM ATI ANO AT3 

CALL MATMUL (ATAPT ,A3T1*AT1»3»3»3) 

CALL MATMUU ATAPT , A 3T3 , AT3, 3,3,3) 

FORM D 

D=2 .0D0*DC0-2 « 000 
FORM DDR 1 
DO 33 1-1,3 

33 DDRU I,U=R3( I , 1 ) -2 . 0D0*R2 1 1 , U *R1 II • 1 1 
FIND SRI 

CALL MATMUL (ATAP#DDR1,D2RI, 3,3,1) 

SRII 3, 11=0. 0D0 
SRII 1, 1I=D2R1( 1,11/D 
SRI ( 2 , II =02R1 ( 2, t ) /D 

FIND RO 

CALL MATMUL ( ATAPT , SR 1 , SR I 1,3,3, 1) 

DO 35 1=1,3 

35 ROI I,1)=R2U ,n-SRIl(I,i) 

FIND V 

CALL MATMUL ( ATI , SRI, SR 11, 3, 3,1) 

00 34 1=1,3 

34 V( I* 1 ) — (— Rl( I,1)*R0{ I, 11 + SR 111 I ,1) )/DT 
FIND SR2 AND SR3 

DO 36 1=1,3 

SRI2I 1,1 )=R2( I,2)-R0( 1,1) 

36 SR I 3( I , 1 ) =R2 ( l , 3)-R0( I » 1 1 

CALL MATMUL! ATAP,SRI2,SR2,3,3,l) 

CALL MATMUL (ATAP»$RI3»$R3,3,3*1) 

150 FORMAT (3(2X,6D20.8«/)) 

PRINT 114 

114 FORMAT! 12X,*W», 17X, *V»,17X, *R0' , 19X , ' SR 1* , 17X, « SR2 • , 17X 
• « * SR3 ( */ I 

PRINT 150* (W( 1,1) ,vn,ll,R0U,n,SRl(!,U,SR2U, 11 ,$R3t 
.1,11*1=1,3) 

DO L 2 L 1 = 1,3 
DV( 1,1 1=999.0 
DRO ( I , 1 ) =999. 0 
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DO 121 J- l *3 
DSR1 l , JI =999.0 

121 ORDQT ( I,J»=999.0 
00 115 1*1,3 

[FIDABSIVSd.lt »*LT.l.0D-12)G0 TO 117 

ovUti>*mi,i>-vsti,in/vs(ifi) 

117 IFtDABStROSJ I . 1> » .LT.1.00-L21GO TO 118 
DRO ( I,l)*IROd, It-ROSI 1,1) t/ROS( 1,1) 

118 IF! DABS( SRRSd , 1 1 ).LT. 1.0D-12)GO TO 119 
DSR< I ,1)=!SR1(I ,i)-SRRS( I ,1) ) /SRRSII *1) 

119 IF(0ABSlSRRSd.2) ) . LT. I .00- 12 ) GO TO 122 
DSR< 1 ,2I = (SR2( I * 1 l-SRRS! 1,2))/ SRRSd .2) 

122 LF I OABS( SRRS( 1,3) ) .LT. 1 • 0D-12)G0 TO 123 
0SRI1 ,3t = CSR3d,l)-SRRSd,3t)/SRRS(I,3) 

123 CONTINUE 

00 115 J-1,3 

IF(OABS(RDOTSd * J I) . LT. 1 . 00-12 ) GO TO 115 
ORDOT ( I , J) = ! RDOTI I , J ) -ROOTS ( I , J ) > /ROOTS d, J ) 

115 CONTINUE 
PRINT 50 

PRINT 110, I DRO < I , 1) , (DSR( I * J) » J = l , 3) ,1*1,3) 

PRINT 51 

PRINT 110, IDVd ,1 ), I ORDOT d,J),J*l,3),I*l,3) 

110 FORMAT I5X»4020«8,/) _ 

50 FORMAT (//,IAX, # DR0*,16X,*DSR 1*,15X,*0SR 2'*15X,'0SR 3 

5 1 * FORMAT I//»14X»'DV*»16Xf * ORDOT l • , 13X, * ORDOT 2S13X, 
.'ORDOT 3',/) 

RETURN 

END 


SUBROUTINE STATE (RR* SRR * SRRI * ROOT ,WI,DPSI»TC,T2,NC*!VIS 
• ,RHO, RHOOtOTHETA. DPHI ,NN* V,RO,SSR ) 

IMPLICIT REALMS ( A-H, 0-Z ) 

DIMENSION ROOT I 3, 3 )»RR< 3, 3) , SRR 0,3) »SRRI I 3,3 1, Nil 3) 
DIMENSION V ( 3 ) , CRO ( 3 ) , APH I ( 3 , 3 ) » ATHETA ( 3, 3 ) , APS II 3 , 3 > , 
1APSIHI 3,3 ) ,A0< 3,3) ,A5(3,3 ) , A6I 3 , 3 ) »CR I 3 ) , SR I 3, 1 ) 
DIMENSION A3T(3,3),R(3,1) ,RA ( 3 , 1 ) , SRP I 3, 1 ) , SRPI 13,1) 
DIMENSION IVI3) ,A7(3»3) ,ATHTAH(3?3),SRI (3,1) » WB I 3 ) 
DIMENSION RHO I 3 t » RHODI 3 I 
DIMENSION D( 3,3 I «Ef 3) 

DIMENSION WBT (3,3) ,WX$RA!3*3) ,RDD (3, 3), RD 2(3, 3) 
DIMENSION R0(3,l) ,SSR(3,3) 

DIMENSION SSRAI 3,3) ,A6T (3,3) 

DIMENSION ABI 3, 1 ) ,8C (3,3) 

DIMENSION XI ( 3) » WO ( 3 ) , HOT ( 3 , 3 ) , WDXSRR 1 3 , 3 ) , SRRA I 3 , 1 ) 
ECU! VALENCE! THETAH,THTAH> 

IV1 1)=5 

I V! 2 ) = 3 
1 V!3)=l 
K3-0 

II = 0 

IF(NC.NE.i) GO TO 6 

CALL DATA! WXO,WVO,WZO, PS 10, THTAH.PS IH,V*CRO, SRR, SA, SB, S 
•C,NN ) 

WYO=. 00005 
THETAH=30.0 
PSI H= 30. 0 

WRITE16,999)SA,SB,SC 

999 FORMAT! 7X, 'SA=' ,F16. 8,7X, * SB= ' , F 16. B , 7X , ' SC= ' ,F16.8,//) 
SAS=SA*SA 
SBS=SB*SB 
SC$=SC*SC 
A= ( SBS+SCS ) 

B=( SAS*SCS) 

C=( SAS+SBS ) 

SIGX = A/C 
SIGY = B/C 
DO 920 1*1,3 
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920 SRR I 3 » I)=SC*DSQRT { 1 • 000-SRR ! l , I > **2 /SAS-SRR ! 2 , ! )**2/SBS 
. ) 

PRINT 47,WX0, WYO, WZO ,PS I 0, THTAH , PS I H 

47 FORMAT (7X, • WXO= * , E 16. 8 , 6X , • WYO= • , E 1 6. 8, 6X, • WZO* • , E 16. 6* 
.//,7X, •PSI0=',616.8,5X, *THTAH=« ,E16. 8,4X, »PSIH* • ,g 16.8, 
*//> 

PRINT 4e,VUI,V<2),V<3) ,CR0tl),CR0(2) ,CR0(3) 

48 FORMAT { 7X, 'XD=* ,E16.8,7X, » YD= • , E 16 . 8 , 7X , * Z0= ' ,E 16. B, // , 
• 7X» 1 X = • ,E16.8,8X, 'Y=*,E16.8,8X, * l- * ,E16.8*//) 

00 910 1=1,3 

910 PRINT 49,I,SftRU*tl,$Rft(2,n,SRR(3,n 

49 F0RMATI7X, 'POINT • , 1 1 , 4X , • SX= • , 616. 8 , 7X , • SV= * , E 16. 8 , 7X , 
.«SZ=* ,£16.8,//I 

PI = 3. 14159265 

PSIO=PSIO*PI/180.O 

THT A H=THTAH* PI/180.0 

PSIH=PSIH*PI/180.0 

W1P = WXO/WZO 

W2P = WYO/WZO 

W3P = 1.0000000000 

PSIP = PSIO 

S 162 = l. 00000000 

TP=.5*(SIGX*W1P**2+SIGY*W2P**2 + W3P*'*2) 
HP=(SIGX**2*W1P6*2+ 5 IGY**2 *W2P**2+ W3P**2 ) +*. 5 
SK=( ( I2.*TP-HP**2I*(SIGY“SIGX) )/< !HP**2-2«*SIGX*TP)*( 1. 
.-SIGY J ) )**.5 
CALL 0CFL1 (RES, SK, IERJ 

XIAMPM (HP**2-2.*S1GX*TPJ * l 1.-SIGYI /($IGX*S I GY I I **.5 
Al=( SIGY*U.-SIGY>/($IGX*( i.-SIGX 11)**. 5 
A2=(SIGX*( l.-SIGX)/! SIGY*( l.-SIGYH )**.5 
A 3= ( ! SIGX*SIGY1**.5*(SIGY-SIGX) ) / < ! l.-SIGX) * { l.-$ IGYH * 
.*.5 

ROW=( (HP**2-2.*SIGX*TP)/I 1 . -S I GX )) **. 5 
DELTA=!$lGY*($IGY-SIGX)/( l.-SIGX) ) * ( W2P**2/R0W**2 I 
WRITE (6,3) Al,A2, A3, ROW, DELTA 
3 FORMAT ( 1X,5(2X,012.5) I 
TAUP=0.0 
T=0 • 0 

DTAU=O.OOOl*RES/XLAMP 
WRITE(6,545)TP,HP,RES,SIGX,SIGY,SIG2 
WRIT6(6,546)SK, XL AMP 

545 FORMAT (6(3X,D12.5) ) 

546 F OR MAT(2(3X, 012*5)*/) 

DO 5 I = 1,3 

DO 5 J = 1,3 
A THE TA ( I , J ) - 0.0000000 
APS 1 1 I, J) = 0.00000 
ATHTAH( I , J ) = 0.0 
APS IN 1 1 , J ) = 0.00000 

5 APH [ < I* J) = 0.00000 

ATHTAH 1 1,1) = DCOSITHTAH) 

ATHT AH( 1 , 3 ) “ -DSIN(THTAH) 

ATHTAH! 2, 2 ) = i. 0000000 
ATHTAH! 3, 3 ) = ATHTAH 1 1,1) 

ATHT AH ( 3, 1 ) “-ATHTAHI 1,3) 

APSIH! 1,1) = DCQS(PSIH) 

APS I H ( 1 , 2 ) = DSIN(PSIH) 

APSIH (2,2) = APS IH( 1, 1 ) 

APS I H( 2, 1 ) =“APS1 H! 1 , 2 ) 

APSIH! 3,3) * 1.00000000 

CALL MATMUL! ATHTAH, APSIH, AO, 3, 3, 3) 

6 IFJT.GT.TOGO TO 901 
IV1S=1 

1 2 = 0 

T = TAUP/WZO 
SCK=1 .-SK**2 
IF! SK.LT. 1.00“05!K3=2 
IF!K3.E0.2)TAUP=TC*WZ0 
IF!K3.EQ.2)T=TC 
■ ARG=XLAMP*TAUP 
CALL DJELF(SN,CN,DN,ARG,SCK) 

WX=( W1P*CN-A1*I W3P/ROW)*W2P*SN*DN)/! 1 • -DELTA*SN**2 ) 

WY=! W2P*CN*DN+A2M W3P/R0WI *W1P*SN) / ( 1 ,-DELT A*SN#«2 I 
WZ=! W3P*DN-A3*W1P*( W2P/R0W)*SN*CN>/( l.-DELTA*SN**2) 
CT=DSQRT<HP*HP-SIGX*SIGX*WX*WX)/HP 
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ST=-SIGX*WX/HP 

SP=SIGY*WY/(HP*CTI 

CP=$IGZ*WZ/(HP*CTI 

SPST=SP«ST 

CPST=CP*ST 

PSIP=PSIP+HP*(SIGY-1. )*SP*SP*DTAU 

PSl=PSIP+HP*TAUP 

TMETA = DATAN2 I ST ,CT I 

PHI = DATAN2 I SP »CP I 

CCT=CT*CP 

SCT=DSQRT{ i.ODO-CCT*CCT> 

CTHETA=DATAN2(SCT,CCT) 

DCTHET=CTHETA*10O. 000/3. 141592 

DPH I =PHI* 180. 00000/3. 14159 

OTHETA = THETA*180. 0000000/3. 14159 

OPSI = PSI*180. 000000/3. 14159 

00 222 I = 1,3 

DO 222 J = 1,3 

APHI I I, J» = 0.000 

ATHETAI I , J ) = 0.000 


IFIK3.LT. 2IGD TO 222 
CT=1 .ODO 
ST-O. 000 
CP= l . ODO 

SP=0 . ODO 

“ 

222 APSI I I ,J) = 

0.000 

APH 1(1,1) = 

1.0000000 

APH 112,2) = 

CP 

APH 1(2,3) = 

SP 

APH I 1 3 , 2 I = 

-SP 

APH I ( 3 , 3 I = 

CP 

ATHETAI i, i) 

= CT 

ATHETAI 1,3) 

= -ST 

ATHETAI 2,2) 

= i. ooooooo 

ATHETAI 3,31 

= CT 

ATHE T A I 3 , 1 1 

= ST 

APSI (1,1) = 

DCOSIPSI 1 

A PS I l 1 , 2 ) = 

DSINCPSI > 

APS I ( 2 , 1 I = 

- A PS I 1 1 , 2 1 

APS I I 2 , 2 ) = 

APSI (1*1) 

APS I I 3 , 3 I = 

1.000 

CALL MATMUL! 

ATHETAfAPSI «A5,3,3,3) 

CALL MATMUL ( APHI ,A5,A6,3,3,3) 

CALL MATMUL 

( A6, AO, A7, 3,3,3) 

00 33 I * 1, 

3 

DO 33 J = 

1.3 


33 A3T( ! , Jl = A7 1 J 1 1 l 
00 949 KK = 1, 3 
DO 950 1=1,3 
950 SRI I,1I = SRRU,KK> 

CALL MATMUL (A3T,SR,SRl,3,3,l) 

DO 10 I * 1,3 

CRin * vi ii*t ♦ croi i > 

10 R 1 1 , 1 ) =CRC II ♦SR I ( 1 , 1 ) 

CALL MATMUL I A7,R,RA, 3, 3, 1 1 

RH01 = 0SQRTIR( 1, l)*R( 1 , 1 ) +RI 2 , 1) *R I 2, ll+RI 3, 1 1 *R I 3 , 111 

RH0IKK)=RH01 

XLA=RAI l,il/RH01 

XMA=RAI2, 1I/RH01 

XNA=RA C 3, 1 l/RHOl 

A4=XLA*XLA/SAS+XMA*XHA/SBS+XNA*XNA/SC$ 

B4=( A4-XLA*XLA/SAS)*SR(1,1 )-XMA*XLA/SBS*SR( 2, 1)-XNA*XLA 
./SCS*SR(3, LI 

SKP I 1. 11 = 2.00000*84/ A4-SRI 1,1) 

SRPI 2,1I = (XMA/XLA)*|SRPI 1 , 1 1 -SR 1 1, 1 1) +SR I 2 , 1 1 
SRP I 3,11 = (XNA/XLA)*(SRPI 1,1)- SR I 1 , 1 1 I + SR I 3 , 1 1 
CALL MATMUL I A 3T , SRP , SRP I , 3 , 3 , 1 1 

RH02=DS0RT( (CRI 1)*SRPI(1, 11 I ** 2 + I CR 1 2 >♦ SRP I (2,il 1**2 + 1C 
•R(3)*SRPI 13,111**2) 

I F I RHD2.LT . RHQ1 I I VI S=I VI S + I VI KK ) 

50 CONTINUE 
948 CONTINUE 

DO 949 I J=l,3 
RRI IJ,KK)=RI I J, II 


19 



949 SRR ! ( I J * KK ) = SR I 1 I J , 1 ) 

WX = WX*WZO 
WY = WY*WZO 
WZ = WZ*WZO 
WB(1)=HX 
W8 ( 2 ) = WY 
WB ( 3 ) = WZ 

xnn = <scs-$BS»/A 

XI(2>=ISAS-SCS)/B 

xn3) = «ses-SAS)/c 

WD( l)=WB( 2>*WB( 3)*XI (1) 

WD( 2)=WB(3)*WBI H*XI ( 2 > 

WD( 3 > -MB ( l >*WB( 2 > *X I 1 3 > 

00 908 J = 1 f 3 
00 908 1=1,3 

908 ROOT ( I ,J)=V( U+A3TJ I , 11 * ( -WB ( 3 J *SRR ( 2 , J ) + WB ( 2 I *SRR < 3 * J ) 
. )+A3TU,2l*<WB< 3)*SRR(1»JI-WBI 1)+SRR(3,J) I + A3T ( I ,3 1 *( -M 
.8<2»*SRR< l,JI*WBm*SRR(2,J)l 

00 899 J = 1 , 3 

099 RHODI J ) = ( RR ( 1 # J ) *RDOT ( 1 , J ) +RR 1 2 , J ) *RDOT ( 2 * J ) *RR ( 3 » J ) *RD 
•0T( 3, J) »/RHO( J) 

DO 909 1=1,3 

909 WHI )=WX*A3T( [,1 J*WY*A3T( I,2)+WZ*A3TU«3> 

930 CONTINUE 

109 TAUP = TAUP + DTAU 
T=TAUP/WZO 

11=1 
T2 = T 

I F ( K 3 • EQ. 2 IGO TO 901 

110 CONTINUE 
GO TO 6 

901 CALL T ILOE ( WB, WBT ) 

IF( T.GT.TC)T=(TAUP-OTAU)/WZO 
CALL MATMUUWBT,SRR,WXSRA,3,3,3) 

CALL MATMUL(WBT,WXSRA, ROD, 3,3,3) 

CALL OOTPROl W8,WB,WW1 
t F ( K 3 . EQ* 2 > GO TO 699 
00 601 1=1,3 
SSRAI 1 , U = SRR( 1,1 ) 

601 R0( t, n=CRU)+A3T( I ,3)*SRR(3,1) 

CALL TILDE (WDyWDT) 

CALL MATMULt WOT y SRR , WDXSRR , 3 , 3, 3 I 
00 606 1=1,3 
00 606 J= L* 3 

606 RDD < I , JJ=RDD(I, J)+WDX$RR< I, Jl 
GO TO 711 

699 CONTINUE 

DO 700 1=1,2 

700 SSRA( 1,1) =-R00 ( I * 1 ) / WW 
SSRA < 3,1) =0.000 

RDD ( 3 » 2 ) = RDD ( 3, 2) -R 00(3,1 I 
ROD (3,31=800(3,3) -ROD (3*1) 

R00(3,l)=0.000 

CALL NATMUL (A3T,R0D,R02,3,3»3) 

DO 703 I = L , 3 

703 R0< 1 , 1 >=RRI 1,1) ♦RD2( I , 1 ) /WW 
711 CONTINUE 

DO 702 J=l , 3 
DO 702 1=1,3 
702 A6T ( I , J ) =A6 ( J , I ) 

DO 705 1=1,3 
DO 705 J = 2 * 3 
705 SSRA ( I,J)=SRR( I « J ) 

SSRA(3,2>=SRR(3,2)“SRR(3, 1) 

SSRA(3»3)=SRR(3,3 ) -SRR ( 3 , l ) 

CALL MATMUL ( A6T , SSRA ,SSR,3,3,3) 

PRINT 112,DPHI,DTHETA,0CTHET 

112 FORMAT (5X, , PHl=*,D16.8r5X, 'THETA 3 ' , 016* B, 5X, * OCAPTHETA* 
• * » D16. 8 , / ) 

911 CONTINUE 
RETURN 
END 
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SUBROUTINE DATA ( Hi , W2, W3 , PS IOI . THTAH1 , P SI Hi , VI . CR 1 , SRRL 
• » SI * S2 * S3 , NN ) 

IMPLICIT REAL*8 ( A— H , 0-2 > 

DIMENSION CRH3)»V1( 3) ,SRR1<3,3 ) 

DIMENSION VP( 3 > ,CR0P(3> * SRRP(3»3> 

I F I NN .GT* l )G0 TO 10 

REAOI 5. iOOJWXOP, WYOP , HZOP , PSI OP , THTAHP , PSIHP * VP ( 1 I • VP (2 
. ) » VP ( 3 I * CROP ( 1) .CR0PC2) .CR0P13) , SRRP (1 • l ) , SRRP l 2 . It t SRRP 
. I I* 2 J , SRRP (2.2) »SRRP( 1,3) ,SRRP< 2,3) 

100 FORMAT ! 5F16* 0 ) 

READ 10l.SAP.SBP.SCP 

101 FORMAT ( 3F 16. 0 ) 

SRRP ( 3,1>=0.0D0 
SRR P I 3 » 2 I =0.000 
SRRP(3,3)=O.ODO 
NN=NN*1 

10 CONTINUE 
S1=SAP 
S2=SBP 
S3=SCP 
PS1HI=PS IHP 
PSI Ol = PS IOP 
THT AH L = THT AHP 
H1=WX0P 
W2=WY0P 
W3=WZ OP 
DO 20 J=l,3 
CRH J l=CROP( J) 

V1IJI =VP ( J ) 

DO 20 1=1,3 

SRR 1 ( I » J I = $RRP ( I » J) 

20 CONTINUE 
RETURN 
END 


SUBROUTINE VISIB(IVIS) 

IF( I V IS.NE* 1 ) GO TO 19 
PRINT 3 

3 F0RMATl/»6X, 'ALL POINTS VISIBLE',/) 

19 IF( IVIS.LT.6) GO TO 17 
PRINT 6 

6 FORMAT!/, 6X, 'POINT NUMBER i NOT VISIBLE',/) 
I VI S= IV I S-5 

17 IFUVIS.LT.A) GO TO 7 
PRINT 8 

8 FORMAT!/, 6X, 'POINT NUMBER 2 NOT VISIBLE',/) 
IVI S= I V I S-3 

7 IF! I VIS.LT.2I GO TO 9 
PRINT 18 

18 FORMAT!/, 6X, 'POINT NUMBER 3 NOT VISIBLE',/) 

9 CONTINUE 
RETURN 
END 


SUBROUTINE DOT PRO ! 0 , E , DOT ) 

IMPLICIT REAL * 8 ! A-H,0-Z ) 

DIMENSION D! 3,1 ) • E ( 3 * 1 ) 

DOT=D!l»l)*E(l,l)+D(2,l)*E(2,l)+D(3»l)*E(3»l) 

RETURN 

END 
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SUBROUTINE MATMUL ( A ,8«C , L »M,N) 

IMPLICIT RE AL*8 ( A-H » 0-Z ) 

DIMENSION A(L,M) ,B(M,N) ,C(L,N) 

DO L 5 I = 1 , L 

DO 15 K= 1 , N 

TEMP=0.0 

DO 14 J= 1 , M 

14 TEMP-=TEMP*A( l » J ) *B ( J , K ) 

15 C( I » K I =TEMP 
RETURN 

END 


SUBROUTINE TILDE ( E, D) 
IMPLICIT REAL*8(A-H,0~Z> 
DIMENSION DI3.3) . E ( 3 ) 

DO 5 1=1,3 
5 0(1,11*0.0 
D( 1 • 2 »*— E (31 
0(2,1 » = E ( 3 ) 

DI 1,3)= E ( 2 ) 

0(3, 1 )=-E(2) 

0 ( 2 , 3 ) =-6 < 1 ) 

D t 3 » 2 ) = E ( 1 1 

RETURN 

END 


Listings of subroutines DJELF and DCEL 1 appear in the listing of program » DUMRKR. 
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DUMRRR Program Listing 


IMPLICIT REAL *8 ( A-H , O-Z I 

DIMENSION 6A ( 11 * 11) , B( 11 • 11 ) ,H( 2, 111 ,R1I3, 1) »OZ (2, l > »R0 
.13*1) , VO (3,1) ,CR( 3, 1) *DX( II *11 .WEIGHT 1 11, U) ,PR0( 11*1), 

• A ( 3 * 3 ) ,C( II, l) , ATI3.3) 

NN=l 

NM=0 

READ! 5, 1 ) TC, DT 

1 FORMAT ( 2F 1 6. 0 ) 

CHECK=l.ODO 

T = TC 
NC= 1 

ROEST -I *0D05 
HST= 1 . 0U-04 

READ < 5»13)PSIH»THETAH,WE*R0(1»1 ) »R0(2,I ) 

READ I 5,13>R0(3,ll .VOU.U »V012,1>,V0<3,I),R1U,1> 

READ! 5, 13 ) R1 (2, l) 

13 FORMAT ( 5F1&. 0) 

Rl(3»l)=0. 0D0 

PRINT 937, PSIH, THETAH, WE ,R0 ( 1 , 1) ,R0 ( 2 , 1) , RO ( 3 , I ) , VO (1 , 1 

. ) , VO ( 2 , L ) , VO ( 3 , 1) , Rl ( 1, I ) , RI( 2, I ) ,R1 ( 3 • 1 ) 

937 FORMAT (5X, • PS I H= • , D16 . 8 , 5X , * TH£TAH=* ,D16.8»5X,*WE-' ,D16 

• .8,/, 5X , * XO- ' ,016.8, 5X, * Y0= « ,01 6. 8 , 5X , ' Z0 = * ,D16.8,/,5X, 
. •XOD=* ,D16.8,5X,« Y00=* , 01 6. 8 , 5 X , • 200= * , DI6 . 8 , / , 5X , * X 1 = • 
.,D16.8»5X,*Y1=*,D16.8,5X,*ZI=',D16.8,//) 

2 CONTINUE 
NM=NM+ l 

I NM=NM/ 10 

P$lH=PSIH*3.141592/180.0 
T HE TAH=THETAH* 3. 141592/180.0 
DO 30 K=l,6 
P SI =WE*T 
PRINT 947, T 

947 F0RMAT15X, 'T=' ,016.8,/) 

948 CONTINUE 

CALL RHOE $ ( 02 * T » NC* P SI H* THETAH » W6,R0 « VO, Rl » CR » A, P S I , ROE 
. ,RUR0D6,NN,NM,INM) 

DZ( 1, 1 )=DZ( 1,1) /ROEST 

DZI2, i)=DZ(2,l)/R06ST/R0eST/WST 

CALL TRP0SE(A,AT*3,3» 

TND=T*WST 

00 38 1-1,3 

RO I I ♦ 1 )=RU(I ,1)/R0EST 

V0( I,1)=V0(I,1)/R0EST/WST 

38 Ri ( I , 1 )=Rl (1,1) /ROEST 
ROE=ROE/ROEST 

k£=WE / WST 

CALL DERI V(H,RO, V0,R1,WE , PS IH, THET AH , ROE , PS I , TNO, AT , A ) 
IF( K . t Q . 6 ) GO TO 33 
DO 31 1=1,2 
DO 31 J= L, 11 
6 ( 2*K-2+ l , J) =H( I , J) 

31 C<2*K-2M,i)=DZU,l) 

T=T+DT 

00 39 1=1,3 

RO I l,l)=RO(l,l)*ROEST 

VO l I , l )=V0 ( 1,1 ) *RQEST*WST 

39 R 1( I ♦ 1 ) =R 1 ( I,1)*R0EST 
hlE=WE*WST 

30 CONTINUE 
33 CONTINUE 

DU 32 J*1,H 

32 6( 11* J)=H( 1, J» 

cm,n=oza,i> 

PRINT 9L9 

919 FORMAT ( 10X, *C MATRIX*,/) 

PRINT 51, t C I 1,1), 1 = 1,11) 
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51 FORMATS 11!5X,D20.8»/ ) / 1 

949 CONTINUE 

DO 4 1*1, li 
DO 4 J*1,1L 

4 BA( l , J ) =B ( l , J ) 

CALL A37I MS( li, BA > 

CALL MATMUL(BA,C,OX»ll»ll»i) 

CF=0.02 

IF(CHECK.LT.i.0D-4)CF=0.i 
lFtCH£CK,LT.l.OD~6)CF=l.O 
OQ 5 1*1,3 

R0( I,1)=R0II ,i)4DX{ I+3,l)*CF 

5 vo( i , n=von , 1 1 +Dxu+6 , u*cf 
R l( if l )=R1< 1,1) +0X110,1 >*CF 
RU2,n=RU2,l)+DX(ll,n*CF 
PSIH=PSIH+OX 1 1, 1) *CF 
THETAH=THETAH*DX< 2, 1 )*CF 
W£=WE+DX!3,1)*CF 
WE=WE*W$T 

DO 40 1=1,3 
ROt I,i) = ROU,l)*ROEST 
40 VOI [»l)=VUl I,ll*ROE$T*WST 
Rl!l,l)=Rl(l«l) *ROEST 
Rl(2,l)=Ri(2,l) *R06ST 
CHECK=O.OUO 
DO 35 1*1,11 

35 CHECK=C ! I , 1 ) *C ! I , 1 ) *-CHECK 
CHECK=DSQRT!CHECK) 

PS1H=PSIH* 180,0/3.141592 
THE TAH=THETA HP 180.0/3.141592 
PRINT 8,PSIH,THETAH,HE,CHECK 

8 FORMA T(5X, *PSIH=* ,D2Q . 8 , 4X, • THETAH= • , 020. 8 »4X , » WE=* ,D2 
.0. 8, 4X, ’CHECK* * ,020.8,//) 

PRINT 9, (ROt I, l I , VOC I ,1) , 1=1,3) 

9 FORMAT! I 5X, ‘RO 1 ,23X, • VO • , / , 3 < 5X ,020. 8 , 5 X, 020. 8, / ) , // ) 
PRINT 10, (Rl{ 1,1), 1=1, 2) 

10 FORMAT! 1 5X , * R 1 * , / , 2 < 5X , D20. 8 , / ) , it I 
PRINT 1 1 , ( DX ( 1,1), 1*1, ID 

11 FORMAT (20X, * D X MATRIX* ,/,4!5X,3( D20 . 8 * 3X ) » // ) ,//) 

950 CONTINUE 

IF!CHECK.LT.1.0D-07)G0 TO ill 
T = TC 
NC=2 
GO TO 2 
111 STOP 
END 


SUBROUTINE DER I V! H, RO, VO , R l , WE , PS I H, THE TAH, ROE , PS I , T, AT 
. ,A) 

IMPLICIT REAL*3 ( A-H, 0-2 ) 

DIMENSION A(3,3),AT(3,3),H(2,ll),RO!3,l),VO{3,l) 
DIMENSION Rl!3, 1) ,0L(3, 1) ,D2(3,1) ,D3!3, l),D4!i,3) 
DIMENSION 05! 1,3) ,06 ( 1,3) ,C!3,1 ) »R IT 11,3) ,R0T!1,3) 
DIMENSION R(3,1),G31(3,1) ,VOTI 1,3) 

DIMENSION DIO 13, 3), WET I 3, 3), Gil! 1,1) 
AA=Rl!l,l)*DC0S!PSl)"Rlt2»l )*OSIN(PS! I 
BB = RU l,l)«DSIN!PSI I +R1 ! 2, 1I*DC0S! PS I > 

00 6 1=1,3 

R(I,1)=R0( I,1)+V0( I , i)*T 

on i ,11=0.0 

D2( I , 11=0.0 
6 D3t I , l)=0.0 
DU 1, l)=1.0 
02 ( 2 , 1 ) = 1 .0 
D3 {3, 11 = 1.0 

CALL TRPQSE(R1»R1T,3,1) 

CALL TRP0SE!R,R0T, 3, 1) 

CALL TRPOSE! V0,V0T,3,1) 

Cll, U=-AA*DCOSITHETAH)*DS1N!PSIH)-BB*DCOSIPSIH) 
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CI2, 1 )=AA*DCOSI THETAH)*DCOS(PSIH>-BB*DSINIPSIH> 
Cl 3,1) =0.0 

CALL MATMULIROT, C,Gll,l, 3, 1) 

HU,n=GiK i,i»/Roe 

Ctl.l l=-AA*DCO$(PSIH)*DS|N I THETAH) 

C I 2 » i )=-AA*D$lN(PSIH)*D$IN( THE T AH ) 

C(3, 1)=-AA*DC0$(THETAH) 

CALL MATMULIROT, C,G11, U3.lt 
H(i,2l=Gllll»l)/R0t 

CI1. l)=-BB*DCOS(PSIH)*DCOSl THETAH )-AA*OS INI PS IH) 
C(2, 1 )=-BB*DSIN|P$IH)*UCOS( THETAH) +AA*OCOSI PSIH) 
Cl 3, l)=-BB*DSIN(THET AH > 

CALL MATMUU ROT , C ,G 1 1 , 1 , 3 , l ) 

HI l*3)=Gll(l, 1 ) *T /ROE 
CALL MATMUL(A*01,G31,3,3, 1) 

CALL MATMUUR1T,G31,GU, 1,3,1) 

HI 1,4J=(R( 1,1 >+Gll( 1,1) > /ROE 
CALL MATMULI A , D2 , G3 1 , 3, 3 , 1) 

CALL MATMUL IR1T ,G31,Gll,i,3,ll 
HI 1 ,5 ) = (R( 2, 1 H-GU(l.l) )/ROE 
CALL MATMULI A, 03, G31, 3,3, 1) 

CALL MATMULIR1T, G31, Gil, 1 ,3,11 
HI l ,6> = (R(3, ll+GU I 1,11 )/KQ6 
HI l,7)=T*Hll,4> 

HI i » 8 ) = T*H I 1 * 5 ) 

H{ 1,9 ) =T*H( 1 , 6 ) 

CALL MATMULI A ,01 ,G31 » 3 , 3 , 1 1 
CALL MATMULIROT ,G31,G1L, 1,3, II 

Hii,iai»icuii,n*Ri<itm /roe 

CALL MATMULI A, D2,G3l ,3,3, 1 ) 

CALL MATMULIROT, G31.G11, 1,3,1) 

Ha,in=iGii(i,i)+Ri(2,n i/roe 

Cl i, 1 )=-AA*DCOS(THETAH)*DSIN(PSIH)-B8*DCOS(PSIH> 
C 12, 1 )=AA*DCOS( THETAH)* DC OS(PSIH)-BB* OS INIPS1H) 
CI3, 1 >=0.0 

CALL MATMULI VOT,C, Gil, 1*3,1) 

HI 2, 1 )=G11 11,1) 

Cl l,l)»8B*DSINIPSIH)*OCO$(THETAH)-AA*DCaSIPSIH> 
Cl 2, M=-8B*DCOS<PSIH)*DCOS(THETAH>-AA*DSIN(PSIH) 
C 13 , 1 ) =0.0 

CALL MATMULIROT, C, Gil, 1*3,1) 

H I 2 , 1 1 - H I 2 , l > ♦ W E *G 1 1 1 1 , 1 1 

Cl 1,1 )=-AA*DCOS I PSIH )*0S INI THETAH) 

Cl 2, l )=-AA*DSIN (PSIH )*DSI N( THE T AH) 

Cl 3, 1 )=-AA*OCO$< THE T AH I 
CALL MATMULI V0T,C,G11, 1,3,1) 

HI 2, 2 ) =G1 1(1,1) 

CD=DCOS(PSI )*WE*R1 12, i l+DSINIPSI )*WE*R1 1 1,1 > 

C(l« 1 )=DCO$(PSIH)*DSIN(THETAH)*CO 
C(2, 1»=DS INI PSIH) *0S INI THE TAHJ*CD 
C( 3, 1 ) =OCOS( THETAH)*CO 
CALL MATMULIROT, C,GLl,l, 3, l) 
H|2,2)=H(2,2)+Gll(l,l) 

Cll*l J*-Rl<2* 1) 

CI2, 1 )=R1I 1* 1 ) 

C ( 3 , 1 ) -0.0 

CALL MATMULI AT, C,G31, 3, 3,1) 

CALL MATMULIROT, G3L, Gil, 1,3,1) 

H12, 31=011(1,1) 

C ( 1, 1 ) =-B6*0C0S I PSIH )*DCO$I THETAH) -AA*DSIN( PSIH ) 
Cl 2, 1 ) =-Bb*DS IN I PSIH )*QCOS( THETAH |+AA*DCOSI PSIH) 
Cl 3, 1 )=-8B*DSIN( THETAH) 

CALL MATMULI VOT.C.Gll, 1,3,1) 

Hi 2 , 3 ) =H( 2 , 3 ) +T*Gll 1 1 » l ) 

C( U l >=-AA*DCOS( PSIH )*DCaS( THETAH I ♦BB*DSINI PSIH) 
CI2, L)=“AA*OS IN I PSIH )*0C0SI THETAH )-BB*DCOSl PSIH) 
Cl 3,1 >=-AA*DSINITHETAH) 

CALL MATMULIROT, C,GU, 1,3,1) 
HI2,3)=HI2.3I*G11I1,1)*T*HE 
CALL TRPOSEI D1 ,04,3, 1 ) 

CALL TRPUSE(D2,D5,3, 1) 

CALL TRP0SEID3, D6 ,3,1) 

DO 12 1=1,3 
DO 12 J= 1 , 3 
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12 WET l l,J)=0.0 
W6T(l,ZI=-W6 
WET ( 2 . 1 1 = WE 

CALL MATMUL ( WE T.RI.G31, 3.3.1) 
CALL MATMUL (AT.G3L.G31.3.3.1I 
CALL MATMUL (DAfG3l»Gll»l»3»l) 

H ( 2 , A 1 =V0 ( 1 » U + CiUtUl) 

CALL MATMULt D5.G31.G11 , 1 . 3, 1) 
H( 2 ,5 1 =VOt 2, l)*Gll( 1 *1 > 

CALL MATMULt D6,G31fGll,l*3, l ) 

H ( 2 • 6 ) =V0 ( 3, n + GlUltll 
CALL MATMULt AT , R l , G31 . 3 » 3 » 1 ) 
CALL MATMUL (DA»G31»GHtl»3»l) 
H(2,7>=RQTtl,l)+Gll(l.l> 

CALL MATMUUD5.G3l»Gll*l* 3.1) 
H(2,8)=R0T(1.2)+G11< 1» 1) 

CALL HATMUL(D6,G31,Gll*i,3,l) 
H(2,9)=R0Ttl,3)*Giit 1,11 
CALL MATMUH WET ,D 1 . G3i » 3 , 3,1) 
CALL MATMUH AT , G3 1 , G 31 * 3 , 3 • 1 1 
CALL MATMUL(R0T»G31,Gll»l»3,l) 
H (2 f 10 I =G1 1 ( 1« l ) 

CALL MATMUL (AT.D1.G31.3.3.1) 
CALL MATMULt V0T.G31. Gllti, 3,1) 
H ( 2 . 10 > — H < 2f 101+GlKhll 
CALL MATMULt WET, 02»G31 t 3. 3, l) 
CALL MATMUL{AT.G31fG3i»3t3.i) 
CALL MATMUUR0T.G31.Gll. 1*3,1) 
H( 2. 1 1 ) =G1 1(1.1) 

CALL MATMULt AT . 02.G31 . 3. 3 , 1 ) 
CALL MATMULt V0T.G31.G11. 1*3.1) 
H(2, 11 l=H(2, 11 1+GLlt l.l) 

RETURN! 

END 


SUBROUTINE RHOE S t DZ * T * NC . PSI H, THETAH , WE .RO, VO . Rl . CR. A , P 
. SI « ROE .RGRODE.NN.NM, INH) 

IMPLICIT REAL*8tA-H»0-Z) 

DIMENSION RR t 3 . 3 ) *SRR(3*3) .SRRI ( 3 . 3 ) ,RDOT t 3 . 3 ) ,WI t3». 

. RHO ( 3 1 • RHODt 3 1 . RO ( 3 . 1 ) , R L ( 3 . I ) » CR t 3 , 1 > • VO 1 3 . 1 ) * 

• WET ( 3«3)*A(3f3).AT(3.3) » G31 1 3. 1 ) 

DIMENSION OZ ( 2 » 1 1 *G32 ( 3 * 1 ) . VI I (3,1) 

DIMENSION R< 3,1 ),V(3,1) 

CALL STATEtRR.SRR.SRRI .RDOT.WI ,DPS I ,T . T2.NC, I VI S.RHO.RH 
.OD.DTHETA.DPHI .DDD.FFF.NN.NM.INMI 
NC = 2 

00 12 1=1.3 
G31 t I ,1 1=0.0 

12 CR< l * 1 )=R0 I I * 1) ♦VOt I . 1) *T 

Atl,il*OCOS(PS! 1 *DCOS( THE TAH )*OCOS ( PSlH)-OSlNtPSI )*DSIN 
. IPSIH) 

AH, 2 »=DCO$( PS I 1 *DCOSt THE TAH 1 *DS I N( PS IH ) «-0S! N t PS I ) *OCOS 
.(PSIH) 

At 1*3 )=-DCOS( PS1 )*D$ INt THETAH) 

A ( 2 * l ) =-0S INt PS I )*UCOS I THETAH) *OCOS ( P S I H I -OCOS ( PSD *DSIN 
. (PSIH) 

A (2, 2 )=-DSlN(PSI ) *DCOS I THETAH ) *DSI N t P$ I HI+DCOS t P$ l )*DCOS 
.(PSIH) 

A (2, 3>=DSIN( PS I ) *0S INt THETAH) 

At3.lt =DSi Nt THETAH) *DCOSt PSIH) 

At 3,2 1=0$ INt THETAH )*DSIN I PSIH) 

At3,3)=0C0StTHETAH) 

CALL TRP0SE(A,AT,3t 3) 

CALL MATMULt AT, Ri »G31 , 3 , 3 , II 
DO 17 1=1,3 

17 Rt I . 1 1 =CR ( I * l I+G3K 1,1) 

CALL OOTPROt R.R .ROES ) 

ROE=DSORT(ROES» 

DO 15 1=1,3 
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DO 15 J = l » 3 
G321 1,11=0.000 
V/ 1 1 ( 1 , l)=0.0D0 
15 WET ( 1 ♦ J ) =0.000 
WETI 1,2>=-WE 
WtTU,l)=WE 

CALL MATMUL(WET,R1,G32,3,3,1) 
CALL MATMUL(AT,G32.V1I ,3,3, l > 
00 18 1*1,3 

is vu, 1 1 =vo 1 1, i) + vim ,i) 

CALL D0TPRQI R » V *R0R0DE ) 

0Z ( l . 1 I =RHO( l )-R0E 

OZ( 2, 1 ) =RHQD ( 1) *RH0( 1 ) -ROR0DE 

RHODE = R0R00E /ROE 

RETURN 

END 


SUBROUTINE DCEL L(RES,AK,IER> 

SUBROUTINE DCEL 1 

PURPOSE 

CALCULATE COMPLETE ELLIPTIC INTEGRALS 
Of THE FIRST KIND 

DESCRIPTION OF PARAMETERS 

RES -RESULT VALUE IN DOUBLE PRECISION 

AK -MODULLUS I INPUT | IN DOUBLE PRECISION 
IER -RESULTANT ERROR CODE WHERE 
I ER=0 NO ERROR 

I hR= l AK NOT IN RANGE “1 TO ♦! 

DOUBLE PRECISION RE S , AK , GEO , AR I , AAR I 

IER=0 
AR l -Z .00 

GE0M0.5D0-AK 1 + 0. 5D0 
GEO=GEO+GEO*AK 
RES=0.500 
IF I GEO I 1,2, A 
L I EK= 1 

2 RE S= 1 . 075 
RETURN 

3 GEO=GEO*AARI 

A GEO*USORT I GEO ) 

GEO=GEO+GEO 
A AR I * AR I 
AR I =ARI +GEO 
R£S=R£S*RES 

IF(GtO/AAR I -0. 999999995D0 >3,5,5 
5 RES*RES/AR I <=6,2831853071 795865D0 
RETURN 
END 


SUBROUTINE DJEL F ( SN , CN , DN , X , SCK ) 
DIMENSION ARIU2I,GE0( 12) 
SUBROUTINE DJELF 


PURPOSE 

COMPUTES THE THREE JACOBIAN ELLIPTIC FUNCTIONS 
SN, CN , DN ■ 


C 


DESCRIPTION OF PARAMETERS 

SN - RESULT VALUE SN < X ) IN DOUBLE PRECISION 

CN - RESULT VALUE CNIXI IN DOUBLE PRECISION 

DN - RESULT VALUE DNIXI IN DOUBLE PRECISION 

X - DOUBLE PRECISION ARGUMENT OF JACOBIAN 
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ELLIPTIC FUNCTIONS 

SCK - SQUARE OF COMPLEMENTARY MOOULUS IN OOUBLE 
PRECISION 

EVALUATION 

CALCULATION IS DONE USING THE PROCESS OF 
THE ARITHMETIC GEOMETRIC MEAN TOGETHER WITH GAUSS 
DESCENDING TRANSFORMATION BEFORE INVERSION OF THE 
INTEGRAL TAKES PLACE. 


DOUBLE PRECISION 5N,CN, DN , X » SCK , AR I »GEO ,CM, Y , A, B, C , D 

TEST MODULUS 

CM=SCK 

Y=X 

1 F ( SCK ) 3 , 1 , 4 

1 0=DEXPU) 

A=l.D0/D 

# B*A^D 

CN=2 *DO/B 
ON=CN 

A= < D- A I / 2 • DO 
SN=A*CN 

DEGENERATE CASE SCK=0 GIVES RESULTS 
CN X = DN X = 1/COSH X 
SN X = TANH X 

2 RETURN 

JACOBIS MODULUS TRANSFORMATION 

3 0=1 .DO- SCK 
CM=- SCK/D 
U=DSQRT l DJ 
Y=U*X 

A A= L.DO 
DN® 1. DO 
DO 6 1=1,12 
L = I 

ARK I ) = A 
CM=DSQRT(CM) 

GEOI I )=CM 
C=( A+CM)*5.D0 

IFIDABS(A-CM)-1.D-9*A)7,7,5 

5 CM=A*CH 

6 A=C 

START BACKWARD RECURSION 

7 Y=C*Y 
SN=DS I N ( Y > 

CN=OC OS I Y ) 

IF (SN) 8, 13,6 

S A=CN/SN 
C = A*C 

00 9 1=1, L 
K*L- 1 + 1 
B=ARl (K) 

A=C*A 

C=UN*C 

DN= (GE0tKl+A)/(8+AI 
9 A=C/b 

A=1.00/DSQRT(C*C+1.DO) 

IFCSNm.il.il 

10 SN=“ A 

GO TO 12 

11 SN= A 

12 CN=C*SN 

13 l F ( SCK ) 1 A , 2 , 2 
1* A=UN 

DN=CN 
CN= A 
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SN=SN/D 

RETURN 

END 


SUBROUTINE A37IMSIN.A) 

IMPLICIT REAL*8IA-H T 0-2) 

250 FURMATI 1H0.20HMATRIX A IS SINGULAR) 

DIMENSION AUMil ,BUll ,C(11I • I X ! MAX tl 1 > , JX JMAX U 1) 
DIMENSION ICHKllII » JCHM 11) ,D{ll,il) 

ISING=0 

MC=N 

DO 132 J = 1 » N 
ICHKI J)=0 
132 JCHK(J)=0 

145 DO 420 I XA * 1 , N 
AMAX=0. ODOO 

DO 160 1 = 1, N 

IF! ICHKI I ) H46, 146, 160 

146 DO 153 J=1 ,N 

IF! JCHKI J» ) 1 47 » 147,158 

147 I F ( AMAX-DA8S I A ( I , J) ) I 150, 158,158 
150 AMAX = DAB$ ( A I I , J ) ) 

1 MAX- 1 
JMAX= J 

158 CONTINUE 
160 CONTINUE 

ICHKI [MAX ) =1 
JCHK I JMAX ) =1 
IXIMAXI I X A ) = I MAX 
J X J MAX ( I X A ) s JMAX 
IF ( AMAX-l.OE-30) 240,290, 290 
240 WRI Tfc I 6, 250 ) 

I S I NC= 1 
STOP 

290 DUM=A(IMAX, JMAX ) 

DO 310 J = 1,N 
IFUCHKI J) )304, 304,310 
304 A{ IMAX, J )=A( IMAX, J ) /OUM 
310 CONTINUE 

DO 380 1=1, N 
I F ( l-IMAX) 320,380,320 
320 DO 340 J=1,N 

IF! JCHK(J)) 330* 330,340 
330 A(I,J)=AII,J)-A(I , JMAX)*A( [MAX, J) 

340 CONTINUE 
380 CONTINUE 
420 CONTINUE 

DU 430 J=1,N 
J T= JX JMAX I J ) 

DO 430 1=1, N 
430 Dt I» J)=A[ I ,JT) 

DO 440 1=1, N 
DO 440 J=1,N 
440 A < I , J ) =0* 0 
DO 450 1 = 1, N 
450 A( I , I > = L » 0 

DO 560 I X A = 1 , N 
DO 452 1=1, N 
452 B { I )=U( I , I XA ) 

I MAX= I X I MAX( I XA ) 

JMAX= JX JMAX I I X A ) 

DO 460 J»1,N 

460 A( I M A X , J )=A( 1 MA X , J )/BI I M AX ) 

DO 558 1=1, N 
I F I I-IMAX>470,558,470 
470 DO 460 J= 1 » N 

480 AI I,J)=A( I , J I -6 ( II*A( IHAX, J) 

558 CONTINUE 
560 CUNT 1 NUc 
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DO 570 IXA=l,N 
1 MAX= I X I MAX ( IXA > 

JMAX=JXJMAXl IXA ) 

I F f I MAX-JMAX) 565* 570,565 

565 DO 566 J = 1.N 
DUM = A I I MAX « J > 

A ( iMAXt J>=A( JMAX, J) 

566 A { JM AX f J ) = DUM 
DO 569 I = I XA ■ N 

I F ( I XIMAX ( I »-JMAXI569,567#569 

567 I XI HAX ( I > = I MAX 

GO to 570 

569 CUNTINUE 

570 CONTINUE 
RETURN 
END 


SUBROUTINE TR POSE ( A , R , N, M I 
IMPLICIT REAL*8(A-H,0-Z> 
DIMENSION A(ll,Rlll 
I R = 0 

DO 10 I = l ♦ N 
I J= I-N 
DO 10 J=1,M 
I J=IJ»N 
1 R= IR*-1 

10 R ( I R I =A ( I J J 
RETURN 
END 


Listings of subroutines STATE , DOTPRO, TILDE, MATMUL and DATA appear in the listing of program DUMRA. 
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